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Preface 


These Lecture Notes have been compiled from the material presented by the 
second author in a lecture series (‘Nachdiplomvorlesung’) at the Department of 
Mathematics of the ETH Zurich during the summer term 2002. Concepts of ‘self¬ 
adaptivity’ in the numerical solution of differential equations are discussed with 
emphasis on Galerkin finite element methods. The key issues are a posteriori er¬ 
ror estimation and automatic mesh adaptation. Besides the traditional approach 
of energy-norm error control, a new duality-based technique, the Dual Weighted 
Residual method (or shortly DWR method) for goal-oriented error estimation is 
discussed in detail. This method aims at economical computation of arbitrary 
quantities of physical interest by properly adapting the computational mesh. This 
is typically required in the design cycles of technical applications. For example, 
the drag coefficient of a body immersed in a viscous flow is computed, then it 
is minimized by varying certain control parameters, and finally the stability of 
the resulting flow is investigated by solving an eigenvalue problem. ‘Goal-oriented’ 
adaptivity is designed to achieve these tasks with minimal cost. 

The basics of the DWR method and various of its applications are described 
in the following survey articles: 

R. Rannacher [114], Error control in finite element computations . In: Proc. 
of Summer School Error Control and Adaptivity in Scientific Computing 
(H. Bulgak and C. Zenger, eds), pp. 247-278. Kluwer Academic Publishers, 
1998. 


M. Braack and R. Rannacher [42], Adaptive finite element methods for low - 
Mach-number flows with chemical reactions. In: 30th Computational Fluid 
Dynamics (H. Deconinck, ed.), Vol. 1999-03 of Lecture Series, von Karman 
Institute for Fluid Dynamics, Brussels, 1999. 


R. Becker and R. Rannacher [31], An optimal control approach to error es¬ 
timation and mesh adaptation in finite element methods . In: Acta Numerica 
2000 (A. Iserles, ed.), pp. 1-101, Cambridge University Press, 2001. 
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Preface 


R. Rannacher [117], Duality techniques for error estimation and mesh adap¬ 
tation in finite element methods . In: Adaptive Finite Elements in Linear and 
Nonlinear Solid and Structural Mechanics (E. Stein, ed.), Vol. 416 of CISM 
Courses and Lectures, Springer, 2002, to appear. 

R. Rannacher and F.-T. Suttmeier [121]. Error estimation and adaptive 
mesh design for FE models in elasto-plasticity. In: Error-Controlled Adap¬ 
tive FEMs in Solid Mechanics (E. Stein, ed.), John Wiley, 2002, to appear. 

Much of the contents of these Lecture Notes is taken from the above articles but 
new material has also been added. At the end of each chapter some exercises are 
posed in order to assist the interested reader in better understanding the presented 
concepts. Solutions and accompanying remarks are given in the Appendix. For 
these practical exercises, sample programs are provided at: 
http://gaia.iwr.uni-heidelberg.de/httpdoc/Research/software.AFEMforDE.html. 



Chapter 1 

Introduction 


We begin with a brief introduction to the philosophy underlying the approach 
to self-adaptivity which will be discussed in these Lecture Notes. Let the goal of 
a simulation be the accurate and efficient computation of the value of a func¬ 
tional J(u) , the ‘target quantity’, with accuracy TOL from the solution u of a 
continuous model by using an approximative discrete model of dimension N : 

A{u) = 0, Ah(uh) = 0. 

The evaluation of the solution by the functional J( ) represents what exactly we 
want to know of a solution. Then, the goal of adaptivity is the ‘optimal’ use of 
computing resources according to either one of the following principles: 

• Minimal work N for prescribed accuracy TOL , 

N —► min, TOL given. 

• Maximal accuracy for prescribed work, 

TOL —> min, N given. 

These goals are traditionally approached by automatic mesh adaptation on the 
basis of local ‘error indicators’ taken from the computed solution, assuming that 
this can indicate local roughness of the ‘continuous’ solution. The main ingredients 
of this process are: 

• rigorous a posteriori error estimates in terms of data and the computed 
solution employing information about the continuous problem; 

• local error indicators extracted from the a posteriori error estimates; 

• automatic mesh adaptation according to certain refinement strategies based 
on the local error indicators. 
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We will demonstrate by examples that the appropriate choice of each of these steps 
is crucial for an economical simulation. Inappropriate realizations which violate 
the characteristic features of the underlying problem may drastically reduce the 
efficiency and accuracy. 

The traditional approach to adaptivity aims at estimating the error with re¬ 
spect to the generic energy norm of the problem, or the global L 2 -norm. However 
this is generally not what applications need. In the following, we will present a col¬ 
lection of examples for such situations, where one is really interested in computing 
locally defined quantities. 


1.1 A first example: computation of drag coefficient 

In order to illustrate the role of adaptivity in the design of a computational mesh, 
we consider a viscous incompressible flow around a cylinder in a channel with a 
narrowed outlet as shown in Figure 1.1. The flow quantities, velocity v and pres¬ 
sure p, are determined by the classical ‘incompressible’ Navier-Stokes equations 


dtv — v Av + v-Vv + Vp = /, V-v = 0. 

The configuration is two-dimensional, with Poisseuille inflow, and Reynolds num¬ 
ber Re = 50, such that the flow is laminar and stationary. The narrowing of the 
outlet causes a so-called corner singularity of the pressure. 



Figure 1.1: Configuration and streamline plot for flow around a cylinder (Re = 50,). 
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The goal is the accurate computation of the corresponding drag coefficient 

J(v,p) := Cdrag := J n T (2i/T-pI)dds , 

of the obstacle, where S is the surface of the body, D its diameter, U the maximal 
inflow velocity, r := |(Vv + Vt; T ) the strain tensor, and d = (0,1)^ the main 
flow direction. 

In order to control the mesh adaptation process in this simulation, one may 
find good reasons to use either of the following heuristic refinement indicators t)k 
on the mesh cells K : 


Vorticity: rjx := hx\\ V xvhWx 


First-order pressure gradient : rjK •= 

Second-order velocity gradient : tjk := h/cllV^VfcH/c- 


Residual-based indicator. rjK := + h 


1/2 

K 


ThWdK + h/c||V*V^||/c> 


: = Z+I/Awh-Vh-Vvfc-Vpfc, 


^|r : 


l[vdnVh-np h ], 

0, 


n l ^ o\i 

if T C r r igid U 

r r 


-O’ 

vdnVh+nph, if T C Tout 




where n is the normal unit vector and [•] the jump across cell interfaces. 


The vorticity as well as the pressure and velocity gradient indicators measure the 
‘smoothness’ of {vh,Ph }» while the heuristical residual-based indicator addition¬ 
ally contains information about local conservation of momentum and mass. As 
competitors, we additionally consider global uniform refinement and refinement 
using a new approach, called DWR method (Dual Weighted Residual method) 
which uses the same residual terms as in the heuristic residual-based indicators 
but multiplied by weights obtained by solving a global ‘dual’ problem. This method 
will be systematically developed in these Lecture Notes. 

The above test case has been designed in order to demonstrate the ability of 
the different refinement indicators to produce meshes on which the main features 
of the flow, such as boundary layers along rigid walls, vortices behind the cylinder, 
and the corner singularity at the outlet, are sufficiently resolved. Figure 1.2 shows 
locally adapted meshes obtained on the basis of the different refinement indicators. 
The results shown in Figure 1.3 demonstrate that in this case the two ad hoc 
indicators involving only the norm of vorticity or pressure gradient are even less 
efficient than simple uniform refinement. This demonstrates that a systematic 
approach to goal-oriented mesh adaptation is needed which not only takes into 
account local properties of the solution but also the global dependence of the 
error in the target quantity on these properties. 
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Introduction 




Figure 1.2: Meshes with 5,000 cells obtained by the vorticity indicator (left), the 
heuristic residual-based indicator (middle), and the new DWR indicator (right). 



Figure 1.3: Error J(e) in the drag coefficient versus number of cells N, for uni¬ 
form refinement, the weighted indicator obtained by the DWR approach, the vor¬ 
ticity indicator, and the heuristic residual-based indicator. 


1.2 The need for ‘goal-oriented’ mesh adaptation 

Let us illustrate the need for ‘goal-oriented’ mesh adaptation by two further ex¬ 
amples of different types of complexity: a 3-d flow problem where only on locally 
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adapted meshes sufficient accuracy is achieved with acceptable costs, and a 2-d 
flow problem in which the complication arises through the strong interaction of 
mass transport and heat transfer. In both situations, the main problem is the 
generation of error estimates which reflect the local and global dependency of 
the error on the locally observed solution properties. In fact, usually mesh adap¬ 
tation in solving a coupled system of equations for a set of physical quantities 
u = (u l h )f =l is based on ‘smoothness’ or ‘residual’ information like 

n n 

VK ■= Or T] K = Y! t u, K\\ R i( U h)\\K- 

i= 1 i=l 

Here, D\u l h stands for certain second-order difference quotients, and Ri(uh ) are 
certain residuals as introduced in the previous example. Both kinds of indicators 
are easily evaluated from the computed solution and are widely used in practice. 
The proper choice of the weights uj 1 k is crucial for the effectivity of the adaptation 
process. They should include both a scaling due to different orders of magnitude 

4 

of the solution components u 1 , as well as the influence of the present cell K on 
the requested quantity of interest. 

A cylinder flow benchmark in 3-D 

We consider the 3-dimensional flow in a channel around a cylinder with square 
cross section as shown in Figure 1.4. The Reynolds number is Re = 20, such 
that the flow is laminar and stationary. This is part of a benchmark suite for the 
computation of viscous, incompressible fluid flow (see Schafer and Turek [124]). 



Figure 1.4: Configuration of 3D flow around a square cylinder in a channel. 
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The goal is again the accurate computation of the drag coefficient 


J(v,p) := c d 


2 


U 2 DH 


J n T (2i/T—pI)dds, 


where D, H are geometrical parameters, U the maximum inflow velocity and 
d = (1,0,0) T the main flow direction. The desired accuracy is TOL ~ 1% which 
turns out to be a rather demanding task even for this simple flow situation. 

In Table 1.1, we compare the efficiency of global uniform refinement against 
that of local refinement on the basis of a heuristic residual-based indicator and 
the DWR method as already used in the first example. The superiority of mesh 
adaptation by the DWR method is clearly seen. 


Table 1.1: Drag results: a) Q^IQx-element with global refinement, b) Q\/Q\- 
element with local refinement by heuristic residual indicator, c) Q\/Q\-element 
with local refinement by DWR method; the first mesh level on which an error 
smaller than 1% is achieved is indicated in boldface; from Braack et al. [40]. 


a) N 

Cd 

b) N 

Cd 

c) N 

Cd 

15,960 

8.2559 

3,696 

12.7888 

3,696 

12.7888 

117,360 

7.9766 

21,512 

8.7117 

8,456 

9.8262 

899,040 

7.8644 

80,864 

7.9505 

15,768 

8.1147 

7,035,840 

7.8193 

182,352 

7.9142 

30,224 

8.1848 

55,666,560 

7.7959 

473,000 

7.8635 

84,832 

7.8282 

— 

■■El 

1,052,000 

7.7971 

162,680 

7.7788 

00 

7.7730 

00 

7.7730 

OO 

7.7730 



VtouSMt 



Figure 1.5: Refined mesh and zoom into the cylinder vicinity obtained by the heuris¬ 
tic residual-based indicator; from Braack et al. [40]. 
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Figure 1.6: Refined mesh and zoom into the cylinder vicinity obtained by the DWR 
method; from Braack et al. [40]. 


Figures 1.5 and 1.6 show meshes which have been obtained by using refinement 
on the basis of the heuristic residual-based indicator and by the DWR method. 
The heuristic indicator fails to properly refine the area behind the cylinder which 
causes its poorer drag approximation. 


A heat-driven cavity benchmark in 2-D 

We consider a 2-dimensional cavity flow. The flow in a square box with side length 
L = 1 (see Figure 1.7) is driven by a temperature difference Oh — 6 C = 720#, 
between the left (‘hot’) and the right (‘cold’) wall, under the action of gravity g 
in the vertical direction. The Rayleigh number is Ra ~ 10 6 making this problem 
computationally demanding. Here, the quantity to be computed is the average 
Nusselt number (mean heat flux) along the cold wall defined by 


J(u) = (Nu) c : = 



Kd n 0 ds, 


where Pr is the Prandtl number and po, #o are certain reference values for viscos¬ 
ity and temperature. The underlying mathematical model is the low-Mach number 
approximation of the stationary compressible Navier-Stokes equations which is ex¬ 
pressed in terms of the set of primitive variables u = {v,p,0} denoting velocity, 
pressure and temperature. In this case, due to the large temperature difference, 
the usual Boussinesq approximation is not sufficient (see Becker and Braack [24]). 
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Figure 1.7: Flow in heat-driven cavity: velocity norm isolines (left) and temperature 
isolines (right); from Becker and Braack [24]- 


The meshes shown in Figure 1.8 indicate that the heuristic residual-based indicator 
induces mesh refinement mainly in those areas where the velocity is dominant while 
the weighted error indicator obtained by the DWR method puts more emphasis on 
the region along the hot boundary where the temperature gradient is dominant. 
The latter, given the quantity we want to compute, seems to be more important 
for capturing the heat transfer through the cavity. This is confirmed by the results 
shown in Table 1.2 which show that on the meshes generated by the properly 
weighted error indicators the accuracy in computing the Nusselt number is better 

by almost one order. 

Table 1.2: Computation of the Nusselt number in the heat-driven cavity by the 
heuristic residual-based indicator (left) and the weighted indicator by the DWR 
method (right); comparable error magnitudes are indicated in boldface; from Becker 

and Braack [24]• 


N 

(Nu) c 

error 

523 

-8.86487 

1.8-10- 1 

945 

-8.71941 

3.3-10-* 

1717 

-8.66898 

Uif 5 " 

5530 

-8.67477 

1.2-10- ,J 

9728 

-8.68364 

3.010 -3 

17319 

-8.68744 

8.5-10 -4 

31466 

-8.68653 

6.910" 5 








N 

(Nu) c 

error 

524 

-9.09552 

4.1KT 1 

945 

-8.67201 

1.5 10-^ 

1708 

-8.49286 

1.910-’ 

3108 

-8.58359 

l.o-io-’ 

5656 

-8.59982 

8.7-10- Ji 

18204 

-8.64775 

3.910-^ 

32676 

-8.66867 

1.8-10-^ 

58678 

-8.67791 

8.7-10- 3 

79292 

-8.67922 

7.4-10- 3 
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Figure 1.8: Sequence of refined meshes for the heat-driven cavity with N « 
500,5500,56000 cells: heuristic residual-based indicator (top row), weighted error 
indicator by the DWR method (bottom row); from Becker and Braack [24]- 


1.3 Further examples of goal-oriented simulation 

In this section, let us present a collection of further examples from different areas 
in which ‘goal-oriented’ adaptivity is needed and has already proven to be superior 
over other more ad-hoc approaches. Some of these examples will be discussed in 
more detail in the course of these Lecture Notes. 

Example 1.1. CARS signal in a flow reactor (Becker et al. [26]): 

J(p,v,0,w) k / wfcrdo. 

J^'cars 

In chemical flow reactors Coherent Antistokes Raman Spectroscopy (CARS) is 
used to produce signals which reflect the distributions of the mole fractions of cer¬ 
tain reaction products Wi along the measurement line Tcars • The mole fraction 
of the species Wi is determined by the balance equation 

c p d t Wi - V(pDiVwi) + pv-Vwi = fi(0,w), 

together with equations for the other physical quantities, pressure p , velocity v , 
temperature 6 , and density p . The final goal is to determine certain reaction 
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velocities in the chemical process from these measurements. The accurate compu¬ 
tation of the corresponding quantities is of crucial importance for this parameter 
estimation process. 

Example 1.2. Mean surface pressure of a body in an inviscid flow (Hartmann [71]): 

J(p,v,e) / pdo. 

Js 

In the absence of viscosity, this quantity relates to the drag coefficient of the body. 
The density p , the momentum pv and the energy e are determined by the Euler 
equations, setting p = (7 — l)(e — \pv 2 ) , 

d t p + V-(pv) = 0 , 

d t {pv) + S7'(pv®v) + Vp = pg , 

c p d t {pe) + c p V-(pev+pv) = h. 

Example 1.3. Boundary mean stress of an elasto-plastic body (Rannacher and 
Suttmeier [119]), 

J(w,cr) := / n T ando, 

Jv D 

over the ‘clamped’ part T d of the boundary. The stress tensor a and displacement 
u are determined by the Lame-Navier equations with a pointwise stress constraint, 

— V-cr = /, o — Ce, e = ^(Vii+Vw T ), \a\ < gq. 

Example 1.4. Local intensity measurement of a seismic signal (Bangerth [14]): 

rT+6 

J(u):= / u(x ohs ,t)dt. 

JT-8 

The displacement u is determined by the acoustic or elastic wave equation 

pd\u — V-(aVw) = 0 . 


In applications the elastic coefficient a has to be determined from measurements 
of this kind at varying points x Q bs resulting in an inverse problem. 

Example 1.5. Observed light emission of a proto-stellar dust cloud (Kanschat [93]): 


J(u) := t 

J ra* 


u{x , 0obs> ^ 0 ) do. 


obs >0 


The radiation intensity u is determined by the radiative transfer equation 


ro'Vu + (k+p)u — B — 



Due to the distance of the light-emitting object, the observer (located on a satel¬ 
lite) measures only the mean value of the intensity integrated over that part of 
the surface of the object seen by the observer, and at some fixed wave length Aq . 
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Example 1.6. Critical eigenvalue in hydrodynamic stability analysis (Heuveline 
and Rannacher [77]): 

J(p,v, A) := X crit . 

The triple {p, v, A} is determined by the eigenvalue problem of the Navier-Stokes 
equations linearized about some stationary base flow v , 

—isAv + v*Vv + v-Vv + Vp = Av, V'V = 0. 

Example 1.7. Cost functional in an optimization problem (Becker et al. [28]): 

J(it,g) := J cost (u,q). 

The state variable v and control variable q are determined by 

Jcost(u, q) —► min, A(u) + Bq = 0. 

The state equation may, for example, consist of the ‘incompressible’ Navier-Stokes 
equations containing some control parameter q (e.g., boundary control) and the 
cost functional may be the drag coefficient which is to be minimized. 


1.4 General concepts of error estimation 

In the following, we will introduce some of the main concepts of the DWR method 
for ‘goal-oriented’ a posteriori error estimation within the framework of linear 
algebra. We will use the same concepts later for differential equations as well and 
use this simple example only to introduce the most important aspects. 


Traditional error estimation 

For regular matrices A, Ah E E nXn , and vectors 6, bh € K n , consider the prob¬ 
lems of finding x, Xh E R n from 

Ax = 6, A h x h = b h , (1.1) 


where h is a parameter indicating the quality of approximation, i.e., Ah —> A and 
bh —► b , as h — ► 0 . In this context, we introduce the notation of the approximation 
error e := x — Xh , the truncation error r := AhX — bh, and the residual p := 
b — Axh • Usually, a priori error analysis is based on the truncation error, and uses 
the identity 

A h e = A h x - A h x h = A h x - b h = r , 

to derive an a priori error bound involving a ‘discrete’ stability constant cs,h : 



< 


CS,h 


1 



( 1 . 2 ) 
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In contrast to that, the a posteriori error analysis uses the relation 


Ae = Ax — Axh = b — Axh = p, 


to derive an a posteriori error bound involving a ‘continuous’ stability constant: 


e 


< cs 






Notice that the a priori error analysis is based on assumptions on the stability 
properties of the ‘discrete’ operator Ah , which may be difficult to establish for 
the particular approximation, while the a posteriori error analysis uses stability 
properties of the unperturbed ‘continuous’ operator A which are often available 
from regularity theory. Further, the truncation error r is not so easily computable 
in practical applications. 


Duality-based error estimation 

In order to avoid the aforementioned drawbacks and to estimate the error also 
with respect to arbitrary moments of the solution, we employ a ‘duality argument’ 
well-known from the error analysis of Galerkin methods. For some given j G M n 
assume that we want to estimate the value of the linear error functional 

J(e) = J(u)-J{x h ) = (e,j). 

For the determination of this error, consider the solution z G M n of the associated 
dual (or adjoint) problem 



This leads us to an identity about the error 

J(e) = (e,j) = ( e,A*z ) = (Ae,z) = (p,z), 


and finally to the ‘weighted’ a posteriori error estimate 


n 


|^(e)| < V|a 


Zi 


2=1 



In this estimate the residuals pi are easily computable but the determination of 
the weights z^ requires the solution of the auxiliary problem (1.4). The gain in 
using these weights is that they tell us about the influence of the ‘local’ residuals 
Pi on the error in the target quantity J(u ). 

We want to extend the concept of duality-based a posteriori error analysis to non¬ 
linear problems. Let differentiable mappings A(-), Ah{-) : M n —> R n and vectors 
6, bh G M n be given and consider the problems 


A(x) = b, A h (x h ) = bh- 


( 1 . 6 ) 
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Using the Jacobi matrix A f (x) and the nonlinear residual p : 
have the following identities, for an arbitrary y E R n : 


A(x h ), we 


(A(x)-A(x h ),y) 


/ (A'(x h +se)e,y)ds 
Jo 


(■ Be,y ), 


with the matrix 


B := B(x,Xh) := / A'(xh+se) ds. 

Jo 


For any given error functional 
corresponding dual problem 


J(-) = let z G M n be the solution of the 

B*z = j, 


where regularity of B* is assumed. We obtain the error identity 

J(e) = (e,j) = (e,B*z) = ( Be,z ) = {A(x)-A{x h ),z) = ( p,z ), 
which leads us to the weighted a posteriori error estimate 




i 


(1.7) 


i= 1 


The use of the above error estimate requires the evaluation of the weights \z%\. 
However, z cannot be computed since it depends on the (unknown) error e 
through the definition of B. To this end, we approximate the matrix B by 


B(x,x h ) « B := B(x h ,x h ) 


f A'(x h ) 
Jo 


ds = A'(x h ), 


and solve the linearized dual problem, i.e. in practice an approximation of that, 

A'hixtfzh = jh■ 


This leads us to the approximate error estimate 


• %/ 

|J(e)| &Tj:= \{p,z h )\ < ^Wi\pi\, 


i •— \Zfi, 


( 1 . 8 ) 




The error by this approximation can be estimated. With z the solution of 

B*z = A'{x h )*z=j, 

we have by definition: 


|J(e)| = \(e,B*z)\ = \(Be,z)\- 

< \((B-B)e,z)\ + \(Be,z)\ 

< \((B-B)e,z)\ + \{p,z)\ 

< \((B-B)e,z)\ + \(p,z-z h )\+fj. 
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Using the estimate, 

\B — B\\ = f {A*(xh) — A f (xh+se)} ds < ^Z/||e||, 

Jo 

with the Lipschitz constant L' of A !, we obtain 

|7(e)| < rj + fZ/||e|| 2 p|| + ||p|| \\z-z h \\. (1.9) 

Hence, assuming that \\z\\ can be controlled, the approximate error estimator fj 
is the dominant error term. However, notice that the derivation of (1.9) is based 
on typical ‘linear algebra’ arguments as it does not observe a possible dimension- 
dependence of norms and constants. In the PDE context, considered later on, 
we will have to argue more carefully in estimating the effect of linearization and 
approximation of the dual problem. 

Overview 

The further contents of these Lecture Notes are as follows: In Chapter 2, we apply 
the principle of duality-based error estimation to the Galerkin approximation of 
ODEs. The following Chapter 3 is one of the core parts of these Lecture Notes. 
There, we develop the DWR method for the Poisson equation as the prototype 
of an elliptic problem. Then, Chapter 4 is devoted to the practical realization 
of error estimation and mesh adaptation. So far, the development of the DWR 
method is largely based on heuristic grounds, though strongly supported by com¬ 
putational experience. In Chapter 5, we discuss some of the central theoretical 
questions related to the justification of this approach but quickly get to the limits 
of theoretical analysis. Chapter 6 is the second core part of these Lecture Notes. 
It presents an abstract version of the DWR method for general nonlinear varia¬ 
tional problems. This abstract approach is used in Chapter 7 for developing an a 
posteriori error analysis for the Galerkin approximation of eigenvalue problems. 
Another application is presented in Chapter 8, where the Galerkin approximation 
of optimization problems with PDE constraints is considered using the classical 
Lagrangian formalism. In Chapter 9, we realize the DWR method for the space- 
time discretization of nonstationary problems, with the heat and wave equations 
as model cases. Chapter 10 deals with the applications of the DWR method to lin¬ 
ear and nonlinear problems from Structural Mechanics and Chapter 11 contains 
various results on the approximation of the incompressible Navier-Stokes equa¬ 
tions as one of the basic models in fluid mechanics. In the last Chapter 12, some 
current developments and open problems are addressed. 



Chapter 2 

An ODE Model Case 


In the following, we consider the realization of the ideas sketched in the Introduc¬ 
tion for the initial value problem of an autonomous ODE system: 

= /(^(t)), tG/:=[0,T], u(0)=u 0 . (2.1) 

We assume that the function/(•) is Lipschitz continuous and that the solution u 
exists on the interval [0, T] . Very often, in applications, the goal in numerically 
approximating the solution is to know the end-time value 

J(u)=u{T). 

This goal may be reached by using the traditional Finite Difference (FD) method 
or the Galerkin Finite Element (FE) method with suitable error control and step- 
size selection strategies. The material of this chapter is mainly taken from Bottcher 
[37], and Bottcher and Rannacher [38]. 


2.1 Finite differences and finite elements 


We start with a brief outline of the main philosophies underlying the Finite Differ¬ 
ence and the Galerkin Finite Element schemes. Both types of methods are defined 
on time grids on the interval I described by 


0 = £ 0 < •. • < t n <... < t N = T, 




n—15 





t 


n —1 • 


The finite difference method 

We exemplarily consider the (implicit) backward Euler scheme : Find U n ~ u n := 
u(t n ) , such that 


Uq = Uq, 


U n — Un- 1 4- k n f(U n ) n > 1. 


( 2 . 2 ) 
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Assuming the step size k n to be sufficiently small (or /(•) to be negative mono¬ 
tone), the discrete solution U n exists for all n = 1 The corresponding 

truncation error 

T n (u) := ^(Un-Un-i) - /(«„), 


is bounded like 



Viewing the error e n = u n — U n as particular solution of the discrete scheme with 
right-hand side r n , the discrete Gronwall lemma yields the usual a priori error 
estimate 


£n 


N 

<K^k 

71=1 


n 


T 


71 


(u) 


K ~ exp 




where Lf(t) is the Lipschitz constant of /(•) along the exact solution u . On this 
basis, we may use the following explicit formula for step-size control: 


TOL 

KT\\t»(U)\\’ 



where the following approximation of the truncation error is used: 


T n (u) = k n T°(u) + 0(kl) « k n T° 



Remark 2.1. The a priori error estimate (2.3) suffers from two short-comings: 

• Usually the growth factor K = K(T,Lf) is unknown as it depends on the 
exact solution. 


• The step size formula requires an estimate for the ‘leading’ truncation er¬ 
ror r^(u) along the exact solution. This, however, has to be generated (for 
example by local /i-extrapolation) from the computed solution U n . 


The Galerkin finite element method 


We consider the approximation by the so-called discontinuous Galerkin dG(0) 
(finite element) method. Here, the space 

s[ 0) (/) := {<p : / -> R d ,v lIn € Po(I n )h 

consisting of piecewise constant functions, is used as ‘trial’ and ‘test’ space. The 

‘Galerkin approximation’ U E S^ is determined by requiring that C/J" := uo 
and 



{U'-f{U)^)dt + {[U] n -u 




€ Sl 0) 




2.1. Finite differences and finite elements 
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Here, we have used the standard notation 



:= lim <£>(*), 

tUn 






n > 


for left- and right-sided limits, and jumps of possibly discontinuous functions. 
Clearly, the continuous solution u also satisfies the variational equation (2.5), 
and for the error e := u — U , there holds the nonlinear Galerkin orthogonality 
relation 

N r 

(e'-f(u)+f(U),iP)dt + ([e] n - 1 ,iP+_ 1 )}= 0 VtfeSjW ( 2 -6) 

n —1 ^ 1 71 

Since the test functions are allowed to be discontinuous, the globally formu¬ 
lated dG(0) method reduces to a time-stepping scheme which, in the present 
‘autonomous’ case, is actually equivalent to the backward Euler scheme for the 
values U n := U~ : 

Uq — Hq, U n U n —i = Aj n> /*(C/ n ), 77. ^ 1. 


Now, the a posteriori error analysis for the end-time error ||ejv|| via duality ar¬ 
gument proceeds as follows. We consider the (backward in time) dual problem 


-z f -B(t)*z = 0, T>t> 0, z(T) = \\e N \\~ l e N , (2.7) 


with the operator 



f x (U+se)ds, 


or in weak formulation, 

N 

£ 

71=1 



n 


(<£>, —z' — B*z) dt 





Taking </? := e, and using integration by parts and Galerkin orthogonality, we 
conclude that, with an arbitrary Z E S^\ 
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For the special choice Z = z , the interval-wise mean value of z, 



we obtain the error representation 


N 

IMI = 

71 = 1 

From this, we infer the following two types of a posteriori error estimates. 
• Local ‘weighted’ a posteriori error estimate: 




N 

< [£/]„_i 

71=1 




with the ‘interpolation constant’ cj 
represent the sensitivity of the error 
Pn — ll^n 1 P-^]n — 11| • 



The weights u n = f 7 ||z'||ds 

Tl 

with respect to the local residuals 


• Global a posteriori error estimate: 


\cn II < c i 



dt 


cics max {k n p n }. 

l<n<N 


( 2 . 10 ) 


The stability constant cs = f f ||^ / ||d^ represents the ‘global’ sensitivity of 
the error ||e^|| with respect to the maximal residual. 

On the basis of the above estimates, we obtain the following ‘implicit’ a posteriori 
step-size control strategies (k n the old and k f n the new step size): 


cjNknPnUJji or cjCskjiPfi 


> TOL 
~TOL 
< \ TOL 



Here, the weights oj n and the stability constant cs can be approximated by the 
following formula, given a numerical approximation Z oi z: 






dt « ||[Z] n ||, 




Ell^lnll- 


Remark 2.2. The a posteriori error estimates (2.9) and (2.10) also have drawbacks 
which should be pointed out: 

• The evaluation of the weights uj n requires the solution of the dual problem 
over the whole time interval [0, T]. 



2.2. Efficiency comparison: FD versus FE method 
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This ‘global’ step size control is ‘implicit’ since it involves the simultaneous 
adaptation of all time intervals {/i,..., in each adaptation cycle. 


Remark 2.3. Duality-based a posteriori error analysis for the Galerkin approxima¬ 
tion of ODEs has also been considered in Estep and French [59], and Estep [58]. 
Here, the error estimates contain global stability constants which are derived by 
analytical arguments. 


2.2 Efficiency comparison: FD versus FE method 


In the following, we want to compare the efficiency of the step-size selection strate¬ 
gies described above for the FD and the FE method. We emphasize that in the 
present situation (‘autonomous’ ODE), both methods are only different ways of 
writing the very same scheme. Nevertheless, we will see that using either the FD 
or the FE approach gives quite different results. Consider the special scalar initial 
value problem 

u'(t) = u(t) 2 , 0 < t < T <1, u( 0) = 1, 

with the singular solution (see Figure 2.1) 



1 


1 -t 


Let the goal of the computation again be the approximation of the end-time value 
u(T). For the different mesh adaptation strategies designed for the backward Euler 
scheme and the dG(0) method, we want to estimate the work which is required 
for computing with accuracy TOL, in dependence on T —» 1. 




Figure 2.1: Singular solution (left) and corresponding dual solution (right) for the 
evaluation of the end-time error. 
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The dG(0) method 


Let the step-size control be based on the global a posteriori error estimate (2.10), 


|ejv| < c/cs max {fe n p n }, 

l<n<N 


Pn •— k n |[f^]n— 1 


with the interpolation constant 
dual a priori bound 


ci = A and the stability constant cs from the 



/ 


dt =: cs 


o 


The criterion for the choice of the local step size is 


cic s k n pn — TOL k n 


TOL 

ClCsPn 


\cn\ < TOL 


N 


Notice that, by construction, 5Zn=i kn 


T . For the explicit estimation of the 


work count, the parameters in this formula are approximated as follows: 
i) Suppose that the step size is already small enough such that 


P 


n 


kZ 1 1 [U] 


n 


n 


k 


-l 

n 


u n -u n -1 «sup 


U 


(l-*n) 


-2 


n 


ii) The stability constant cs is essentially determined by the dual problem (after 
linearization about the continuous solution u ) 


z'(t) 


2 


l-< 


z(t), 1 > T > t > 0, z(T) = 1, 


with the solution 


z(t) = exp I 2 


r ds 

Jt 1-8 


z(T) 


1 -t 
l—T 


Hence, 





dt < 


1 


(1 -T) 


This leads us to the step-size distribution 


k 


n 




(1—t„) 2 (l—T) 2 TOL 


The corresponding work measure N is determined as 


N 


N 




1 


N 


71= l 


(l-T)2TOL^ 


^ ^ k n {\ t n ) 


-2 


1 




(1—T) 2 TOL 



dt 


(l—T ) 2 ’ 
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and consequently 


N 


1 


1 


(1 — T) 3 TOL 


For later purposes, we note that the local weights are determined by 


UJ 


n 



z'(t)|dt 


k 


1 -t 


n 


n 


n 


(1 -T) 


( 2 . 11 ) 


The backward Euler scheme 

Let the step-size control be based on the global a priori error estimate 


N 


|e/v| < K Yk n \T n (u)\, 


n= 1 


T n (tt)| <\k n SUp 

/ ' 


u 


n 


n 


The corresponding criterion for the control of the local step sizes is 


k 


TOL 


n 


KT sup j | u 


n 


|eiv| < TOL 


The parameters in this formula are approximated as follows: 
i) The Lipschitz constant of /(•) along the solution u is 


Lf(t) 


2 


1 -*’ 


such that the growth factor becomes 



K(T,Lf) « exp [ j Lf(t)dt j « exp 

o / wo 



2 


1 -t 


dt 


(1 -T) 


-2 


ii) The leading truncation error 


T, 


0 


71 


behaves like 


T. 


0 


n 


< \ sup 


u 


It 


(1-tn) 


-3 


n 


This results in the step-size distribution 


fe 


TOL 


n 


TK(T) sup/ I u" 


n 


(1 — T) 2 (l - t n ) 3 TOL. 


The corresponding work count is 


N 


N 




i 


i N 

^vT E 


t r 3 

L n) i 


n=l 


(1—T) 2 TOL ^ 

V 7 71 = 1 


and consequently (compare this to (2.11)), 


N 


1 


1 


(1 —T) 4 TOL 


( 2 . 12 ) 
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Remark 2.4. We note that for both FD-based and FE-based schemes the efficiency 
of the computation can be further improved by using a more localized adaptation 
strategy. This will be considered in the exercises below. 

Remark 2.5. The critical drawback of the heuristic step-size control strategy used 
for the backward Euler scheme is the possible crude under-estimation of the growth 
factor K(T,Lf). On the other hand, estimating K(T,Lf) by an a priori anal¬ 
ysis is oriented at the worst case scenario and usually leads to over-estimation 
rendering the resulting error bound useless. In general, this prevents the step-size 
control strategy based on the a priori error estimate (2.3) from providing a real 
control of the error. In order to overcome this limitation the following two different 
approaches may be used: 

i) The first one is based on the relation 

e n — c n —i 4" k n f {JJn)^n 4" k n T r i(tt) 4" k n O(e n ), (2.13) 

for the error e n = u n — U n , with initial value eo = 0. Using a guess for the 
truncation error r n (U n ) & r n (u ), obtained for example by local extrapolation, 
the solution E n of the linearized error equation 

E n = E n —i 4- k n f f (U n )E n 4- k n T n (U n ), 0 < n < N, (2-14) 

may be used to get a guess for the true end-time error, En « e# • However, this 
does not provide criteria for a time-step adaptation to reduce the error below the 
prescribed tolerance, ||e^|| < TOL . 

ii) An alternative approach employs a duality argument similar to that one used 
for the dG(0) method, but now on the discrete level. Let Z n be the solution of 
the linearized backward-in-time scheme 


Zn—l — Z n 4" k n f (t/n)^n-l, 0 ^ t n <C tjsj. 


(2.15) 


with starting value Zn • Then, using the error relation (2.13), we obtain 

(^n? Z n ) — (^ni Z n Z n _ 1 ) 4" (c n e n — \ , Z n — 1 ) 4" (^n— 1 1 Z u — \ ) 

= -k n (e n , f'(U n )Z n - 1 ) 4- k n (f'(U n )e n , Z n -\) 

“l - ^71(7*71(^0 4 " 0{e n ), Zn— 1) 4 " (^ti—i? Z n — 1)? 


and summing over 1 < n < N, 


N 


(e N ,Z N ) = (e 0 ,Z 0 ) + V] k n (T n (u) +<D(el),Z n - x ) 


n= 1 


For Z N := e^v lle^v 


-l 


and eo = 0, we arrive at 


N 


W^nW < c S,k Vfc„{||r„W|| + 0(|| e „f)}, 


71= l 


(2.16) 


(2.17) 
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with the discrete stability constant cs,k •— maxo< n <iv-i 
quadratic error term, we obtain the error estimate 



. Neglecting the 


e N 




with an approximation r n (U n ) « r n (u) for the truncation error. Here, the gener¬ 
ally too pessimistic a priori bound K{t^^Lf) is replaced by the a posteriori sta¬ 
bility constant cs,k • In practice, cs k will be determined from the computed dual 
solution Z n using again a suitable guess for the starting value Zn — ^n\\^n\\~ 1 • 
The step-size selection strategy based on the estimate (2.18), is generically implicit 
like that proposed for the dG(0) method, but is capable of producing useful error 
bounds. Surprisingly, it has not found much attention in the FD community. 


2.3 Exercises 


Exercise 2.1. Consider the autonomous initial value problem 

u'(t) = u(t) 2 , 0 < t < T <1, u( 0) = 1, 

with the solution u(t) = (1 —t )~ l . Compute u(T) by the backward Euler scheme 

U n = U n —i + k n U 2 , n > 1, Uq — 1. 

Let the step sizes k n be chosen on the basis of the a priori error bound (2.3), 


N 

|ejv| < K(T) Y2 k n T n (u), r n (u) = A;„r°(u) + 0(kl), 

n= 1 

according to the implicit control 

, ( T0L 

kn ~ \nK(T)\t°{u)\ 

What is the asymptotic work count N = iV(TOL, T) as T —> 1 ? 

Exercise 2.2. Consider the model problem of Exercise 2.1. Let in the dG(0) method 
the step-size selection be based on the weighted error estimator (2.9), 



N 


&N ^ C/ ^ ^ k n p n Ld n , 


n =1 


according to the implicit control 



TOL 
N Priori 


What is the work count N = iV(TOL,T) as T —* 1 ? 
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Exercise 2.3. The dG(O) method is only of first order. By the same principle, im¬ 
plicit dG(r) methods of any order r > 1 can be designed. Alternatives are the so- 
called cG(r) (continuous Galerkin) methods which are Petrov-Galerkin methods. 
The simplest representative is the second-order cG(l) method which uses (contin¬ 
uous) piecewise linear trial and (discontinuous) piecewise constant test functions: 


U( 0) = wo, 



e S< 0 ) (I). 


This method is closely related to the implicit midpoint rule 


Uq — Uq, U n — U n —i H - & n /( 2 1 4“ U n }), 77. ^ 1. 


Develop a residual-based a posteriori error estimate for the end-time error 
of the cG(l) method. 




Exercise 2.4 (Practical exercise). Verify the theoretical predictions in Exercises 2.1 
and 2.2 by a computational experiment. Use the backward Euler scheme and the 
dG(0) method for approximating the solution value u(T) in the model problem 

u'(t) = u(t ) 2 , 0 < t < T < 1, it(0) = 1, 

with the solution u(t) = (1 —1) _1 . 

a) Use all the step-size selection strategies developed in the text and in Exercises 
2.1 and 2.2 with the formulas for K(T ), \r®(u) \, p n and uj n as given in the text. 
Determine experimentally the number of resulting time steps N for a decaying 

sequence of tolerances TOL* = 2 -1 , i = 1,2,3,_Monitor the true error e(T) = 

u(T) — Un and compare it with the given tolerance. Interpret the observed results. 

b) Approximate the leading truncation errors |r^| by second-order difference quo¬ 
tients of the computed solution U n , and the residuals p n = , as defined 

in the text. Repeat the above test using these approximations and report the 
differences to the results observed in (a). 

c) Do the same test as in (a) for a sequence of uniform step size distributions with 
k = T/N . You will observe that the performance of equidistant time meshes is 
almost as good as that of the best adaptive procedure based on the ‘weighted’ a 
posteriori error estimate for the dG(0) method. This surprising phenomenon can 
be explained by showing that in this case 

N « (1-T) _2 |log(l-T)|TOL -1 . 

Try to detect this logarithmic behavior in the numerical experiment. So what is 
the benefit of sophisticated local time-step adaptation? 

d) Perform the same experiment with the cG(l) method described in Exercise 
2.3. First, use uniform step sizes and, then, try to use your a posteriori error 
indicator derived in Exercise 2.3. This should demonstrate the superiority of the 
second-order cG(l) over the only first-order dG(0) method. 



Chapter 3 

A PDE Model Case 


In this chapter, we will develop the basics of the DWR method for linear elliptic 
partial differential equations as originally described in Becker and Rannacher [30]. 
As a model configuration, we consider the Poisson equation on a polygonal or 
polyhedral domain ficM d , with Dirichlet boundary conditions: 

-A u = f in O, u\ d n = 0. (3.1) 

The discretization will be by a Galerkin finite element method which is based on 
the variational formulation of (3.1): Find u £ V := Hq(Q.) , such that 

a(u,tp) := (Vu, VVO = (/,V0 € V- (3.2) 


In this chapter, we will exemplarily consider the a posteriori control of the resulting 
approximation Uh with respect to the following goals: 

• Computation of an overview of the solution’s structure (global norm error): 


|V(w— Uh)\\ < TOL, 


u-u h 


< TOL . 


• Computation of ‘displacement’ or ‘stress’ components at some point a G ft 
(point-value error): 


J(u) := u(a), J(u) := diu(a). 

• Computation of ‘mean normal flux’: 

J(u) := / d n uds. 

JdQ 

However, we will always have in mind the accurate computation with respect to 
arbitrary functional values J(u) of the solution. 
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Here and below, we use the following notation: For a domain HcM d , L 2 (Q) 
is the Lebesgue space of square-integrable functions on Q , which is a Hilbert space 
with the scalar product and norm 


(v,w) n 



vwdx, 


v 


n 


<n 



v 


<n 


9 x 1 / 2 

2 dxj 


Analogously, L 2 (dQ) is the space of square-integrable functions defined on the 
boundary dQ equipped with the scalar product and norm 



v 


on 



The Sobolev spaces H l (Q) and H 2 (Q) consist of those functions v E L 2 (Q) 
which possess first- and second-order (distributional) derivatives Vv E L 2 (Ct) d 
and V 2 v E L 2 (0) dxd , respectively. For functions in these spaces, we use the 
semi-norms 



This notation can be extended to Sobolev spaces H P (Q) of arbitrary order p > 1. 
The space H 1 ^) can be embedded in the space L 2 (dQ) , such that for each 
v E there exists a trace v\qsi E L 2 (dCl). Further, the functions in the 

subspace Hq(Q) C H l (Q) are characterized by the property v = 0. By the 
Poincare inequality, 


v 


n<c||Vv||n, v e 


(3.3) 


the H l - semi-norm ||Vv||n is a norm on the subspace Hq(Q.) . If the set Q is 
identical with the domain on which the differential equation is posed, we usually 
omit the subscript Q, in the notation of norms and scalar products, for instance 
v|| = || v ||q . All the above notation will be synonymously used also for vector- or 
matrix-valued functions v : Q, —» R d or M. dxd . 


3.1 Finite element approximation 

The discretization of the model problem (3.1) seeks an approximations Uh E Vh , 
the so-called Ritz projection of u , in a certain finite dimensional subspace Vh C V, 

a(uh,i>h) = (f,iph) Viph&Vh- (3-4) 

The main feature of the Galerkin method for linear problems is the so-called 
Galerkin orthogonality for the error e u — Uh, 


a{e,i>h)=Q, ^h&Vh- 


( 3 . 5 ) 
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The subspaces (finite element spaces) considered have the form 

V h = {v € V : v\K € P(K), K € T„}, 

defined on decompositions of Q into cells K (triangles or quadrilaterals in 
R 2 , and tetrahedra or hexahedra in K 3 ) of width hx — diam(if); we write 
h = max^Gi \ h k f° r th e 9 lobal mesh width. Here, P(K) denotes a suitable 
space of polynomial-like functions defined on the cell K E . In the numerical 
results discussed below, we have mostly used ‘bilinear’ or ‘trilinear’ finite ele¬ 
ments on quadrilateral or hexahedral meshes, respectively, in which case P(K) = 
Qi(K) consists of shape functions obtained via a bilinear transformation from 

A 

the space of ‘bilinears’ Qi(K) = span{l,aq,X 2 ,# 1 X 2 } or ‘trilinears’ Qi(K) = 
span{l,xi,X 2 ,X 3 ,xiX 2 ,X 2 ^ 3 ,X 3 Xi,xiX 2 X 3 } on the reference cell K — [0, l] rf . Lo¬ 
cal mesh refinement or coarsening is realized by using hanging nodes in such a 
way that global conformity is preserved, that is 14 C V (see also Section 4.2). 
For technical details of finite element spaces, the reader may consult the standard 
literature, for instance Ciarlet [46], Johnson [83] or Brenner and Scott [43], and 
especially Carey and Oden [44] for the treatment of hanging nodes. 

We consider the control of the error with respect to some ‘output functional’ 
J(-), i.e., we want to have estimates for the difference J(e) = J(u) — J(uh ). For 
simplicity, we assume here J(-) to be linear. Following the general concept of the 
DWR method, let z € V be the solution of the associated dual problem 

a{(p, z) = J(ip) V</? € V, (3.6) 

and Zh € Vh its finite element approximations defined by 

a(<Ph, z h) = J(<Ph) V<Ph.€V h . (3.7) 

Using this construction together with Galerkin orthogonality, we obtain 

J(e) = a{e,z) = a{e,z—iph) 

= (f,z-iph) - a{u h ,z-ip h ) =: p{u h )(z-tp h ), 4 € 

The so-called residual p(uh ){•) of the Galerkin approximation Uh may be viewed 
as a functional on the solution space V . Cell-wise integration by parts implies 

p(u h )(z-ip h ) = Y2 {(f+^ u h,z-’4’h)K-{d n Uh,z—il> h )dn :} 

Ke Tfc 

= X! {(/+Awh,z-Vvi)/r + \{[dnU h ],z-i>h)dK\m}, 

K€ T h 

where [d n Uh\ denotes the jump of d n Uh across the inter-element edges (in 2-D) 
or faces (in 3-D), i.e., for two neighboring cells K, K f E with common edge T 
and normal unit vector n pointing from K to K f , we set 

[d n Uh\ = [V^ • n] := (Vu h \K'n r — V u h \Knr ) • 
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Since n = — m! , this actually defines the jump of d n Uh across the edge T. For 
later use, we define the cell and edge residuals Rh and , respectively, by 

Rh\K /4-Al/ft, 

r h\r •= 

We collect the previous results in the following Proposition. 

Proposition 3.1. For the finite element approximation (3.4) of the Poisson prob¬ 
lem, we have the a posteriori error representation 

j{e) = ^2 {{ R h,z-i>h)K + (rh,z-iph)dK}, (3.8) 

KeT h 

with an arbitrary 'iph € Vh , and as a consequence the a posteriori error estimate 

|</(e)| < Vu, ■= ^ (3-9) 

K€T/» 

where the cell residuals (‘smoothness indicators ’) pK and weights (‘influence fac¬ 
tors’) ujk are given by 


/ ±[d n u h ], itT<zdK\dCl, 

1 o, if r c dn. 


Pk (II^/iIIk + h K l \\ r h\\dK) ! > 

WK := (Wz-iphWK + h K \\z-xp h \\l K ) 1/2 , 


for an arbitrary iph tVh • 

Remark 3.2. The dual solution z has the features of a ‘generalized’ Green function 
G(K, K') , as it describes the dependence of the target error quantity J(e) , which 
may be concentrated at some cell K , on local properties of the data, i.e. in this 
case the residuals px* on cells K f ; see Figure 3.1. 



? 



p (cell 



e R (cell error) G(K,K') 


residual) 


Green function 


Figure 3.1: Finite element mesh and scheme of error propagation. 
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In order to evaluate the a posteriori error representation (3.8) or the resulting 
a posteriori error estimate (3.9), we need information about the ‘continuous’ dual 
solution z . Since in practice, z is not explicitly known, such information has to 
be obtained either through a priori analysis in form of bounds for z in certain 
Sobolev norms or through computation by solving the dual problem numerically. 
In the following, we will present some examples in which z can be bounded a 
priori or can even be explicitly determined. In later sections, dealing with real-life 
problems, we will always have to rely on the computational approximation of z. 


3.2 Global a posteriori error estimates 


In the following, we will demonstrate how the previous approach can be used for 
deriving the known a posteriori error estimates with respect to ‘global’ norms. 


Example 3.1. ( Energy-norm error) For deriving the usual error estimate with 
respect to the natural energy norm associated with problem (3.2), we choose the 
error functional 

JM = (V^Ve)||Ve||- 1 , 

considering the error e as a fixed quantity. The corresponding dual solution z e V 
is determined by 

o(<p,z) = (V<£, Ve)||Ve|| -1 Vip € V, 

and admits the trivial a priori bound || Vz\\ < 1 =: cs (stability constant). Clearly, 
in this particular case the dual solution is just given by z = e||Ve|| -1 . From 
Proposition 3.1 we infer the estimate 

J(e) = || Ve|| < ^2 Pk^k < [^2 h2 K p K ) 7 ( /l ^ 2 ^) 7 • 

KeT h Ke T h KeT h 


We recall the interpolation error estimate (see, e.g., Scott and Zhang [125]) 


inf ( V {/i 

\l> h ev h \ 1 

Ker h 

to estimate further, 


-2 

K 


z ~^h\\ K + h 


-1 

K 


iph\\ 


dK 



1/2 


< cj ||Vz||, 


||Ve|| < c/^ ^ l|Vs|| < ^ 


1/2 


KeTh 


KeT h 


This results in the classical energy-norm error estimate 


Ve|| < r) E :=c/( h KPx) 


1/2 


KeTh 


(3.10) 


(3.11) 


with pK as defined in Proposition 3.1. 
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Remark 3.3. By Galerkin orthogonality, there holds 

||Ve|| 2 = ||Vit|| 2 - 2(V«, V« h ) + ||Vu h || 2 

= ||Vu|| 2 - 2(V« fc , Vu h ) + ||Vu fc || 2 = ||Vw|| 2 - ||Vu h || 2 , 

such that energy-norm error control turns out to be equivalent to energy error 
control However, this requires the energy form to be a scalar product. 


Example 3.2. ( L 2 -norm error) To derive an estimate with respect to the L 2 norm, 
we choose the error functional 


J{v) = (< P,e) 


-1 


Suppose that the (polygonal or polyhedral) domain Q, is convex. Then, the corre¬ 
sponding dual solution z G VC\H 2 (Q) admits the a priori bound ||V 2 z|| < 1 =: cs 
(stability constant). From the result of Proposition 3.1, we infer the estimate 


< 


PK^K < ' ( Y h K^K) 


1/2 


Kei h 


Kei h 


Ker h 


Using the interpolation error estimate (see, e.g., Brenner and Scott [43]) 


inf f V (A 


-4 

K 


z ~^h\\](+hJ 


z -i>h\\i K }) ' < cr\\V 2 z 


we obtain 


SC,(E 



< 


Ke T h 


Cl ° S { Y h K p K ) 


1/2 


KeT h 


This results in the well-known L 2 -norm error estimate 


< Vl 2 := C/C S ( y; 


Ke T h 



1/2 


(3.12) 


(3.13) 


with px again as defined in Proposition 3.1. In comparison to the energy-norm er¬ 
ror estimate (3.11) the L 2 -norm estimate involves the weighting h^ K which reflects 
its higher order of convergence. 


3.3 A posteriori error estimates for output functionals 

Next, we turn to estimating the error with respect to local error functionals. Let 
TOL denote the accuracy we want to achieve. 
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Example 3.3. ( Point-value error) To estimate the error at some point a G fi, we 
use the regularized functional 


J(u) := \B e \ 1 / udx = u(a) -f 0(e 2 ), 

Jb £ 


where B e is the e-ball around the point a and e TOL. The corresponding 
dual solution z behaves like a regularized Green function, i.e. in 2-D: 


z{x) = g%{x) « log(r(®)), 



x— a| 2 +e 2 . 


By choosing in the estimate (3.9) suitably, the weights have here the form 


LdK 



h 2 K \K\'/ 2 r- K 2 , 



maxrix 

xeK v 


such that 



Vu := Cl y 

K€ T h 




Example 3.4. ( Point-value derivative error) To estimate the error in the derivative 
in direction Xi at some interior point a G Cl , we use the regularized output 
functional 

J(u) \B € \~ 1 f diudx = d{u(a) 4- 0(e 2 ), 

Jb £ 

where again e := TOL. In this case the dual solution behaves like a regularized 
derivative Green function, i.e. in 2-D, 


z{x) = di 9 *{x) 


{x- a)j 
r(x) 2 ’ 


r(x) 


V 


x—a 



and the corresponding weights like 


LdK 



h 2 K \K\ 1/2 r- K \ 


This results in the a posteriori error estimate 


die(a) 



Ker h 




Compared with (3.14), this localizes the region of influence towards the point a 
even more. 
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Figure 3.2: Examples of computed dual solutions: regularized Green function and 
derivative Green function (scaled differently), for evaluating u(a) and d\u(a ). 


Example 3.5. (Mean normal-flux error) Another type of error functionals involves 
integrals over lower-dimensional manifolds. As an example, we consider the com¬ 
putation of the mean normal flux across the boundary, 

J{u) = / d n uds , 

Jon 

where for simplicity is assumed to be the unit circle. The question is: What is an 
‘efficient’ mesh-size distribution for computing J(u)7 Notice that in this simple 
context the computational goal is trivial since it can be reduced to evaluating 
data, J(u) = f Q A u dx = — f Q f dx. However, in more complex situations error 
functionals of this type cannot be so easily computed and are of high interest (c.f. 
the drag coefficient or the average Nusselt number mentioned in the Introduction). 
Here, the corresponding dual problem 

a(<p, z) = (1, d„tp)da V^eVnC 1 ^) 

has a measure solution of the type z = —1 in Q, z = 0 on dQ. Hence, to avoid 
dealing with measures, we use the regularized output functional 

J e (<p) = £ -1 / d n ydx= / d n ipds + O(e), 

Js e Jdtt 

where S £ = {x E Vt : dist{:r, <9f£} < e} and e := TOL . The corresponding dual 
solution is explicitly given by 


— 1 inf l\S £ , 

— £ _1 dist{x, <9f£} in S £ . 
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Figure 3.3: Refined mesh and computed dual solution for the mean normal flux. 


On cells K C Cl\S £ there holds z — IhZ = 0, which leads us to the error estimate 


Js{e) < r?u, : 


E 


Pk^k 


KeJ h ,Kns £ ^(d 


The conclusion is: There is no contribution to the error from cells in the interior 
of Cl. Hence, whatever right-hand side /, the optimal strategy is to refine the 
elements adjacent to the boundary and to leave the others unchanged. In prac¬ 
tice, due to hanging nodes, this may also lead to some refinement in the interior, 
however. 

Example 3.6. [I?-norm error ) Finally, we consider again the L 2 -norm error esti¬ 
mate but for a problem with strongly varying diffusion coefficient, 


-V-{aVu} = f in 0, u\qq = 0. 



This example is intended to show that even when estimating the error in global 
norms, it can be beneficial to keep the dual weights within the error estimator 
rather than condensing them into just one global stability constant. In the con¬ 
sidered situation the dual problem reads in strong form 

—V-{aVz} = e||e|| -1 in H, z\qq — 0, (3-17) 


and the local residuals take the form 


= 5 n ' 

By the same argument which led us to Proposition 3.1 and then to the standard 
L 2 -error estimate, we infer the following two types of a posteriori error estimates: 


Rh\K — f + V-jaV^}, r h i r 































34 


Chapter 3. A PDE Model Case 


• ‘Weighted’ error estimate: 

|e|| < rfli := c/ ^ h 2 K p K UK, uk '■= ||V 2 .z||k- (3.18) 

Kei h 

• ‘Global’ error estimate: 

1 /2 

\e\\ <Vl*'■= cic s y^2 h 4 K px'j , c s := ||V 2 z||n. (3.19) 

T„ 

The residual terms px are defined as before. In both cases, the stability terms 
can be evaluated by replacing the true dual solution by its Galerkin approxima¬ 
tion ||V 2 2 ||k ^ 11V 2 2^11 k , with some second-order difference operator V 2 . The 
interpolation constant is typically of size c/ « 0.2. The error-dependent func¬ 
tional J(-) = (-,e)||e||~ 1 is evaluated by replacing the unknown solution u by a 

patch-wise higher-order interpolation I { 2h^h Of U h , 

r(2) 

e « I± h J u h -u h . 

This gives us approximate L 2 -error estimators denoted by fj^ 2 {uh) and r/L 2 (^/i), 
respectively. We want to compare the performance of these two L 2 -error esti¬ 
mators by a numerical experiment. To this end, consider the particular setting 
Q = (—1, l) 2 and a(x) = 0.1 4- e 3 ( Xl+X2 ^ , with a sinusoidal solution u(x) and 
corresponding right-hand side /. In this calculation the mesh adaptation tries to 
equilibrate the local ‘error indicators’ r\x — h 2 K px^K and r\x = h'xPx > respec¬ 
tively. (This and alternative strategies for mesh adaptation will be discussed in 
more detail in the next chapter.) 



Figure 3.4: Point-value errors obtained by ry^ 2 (left) and rff 2 (right, scaled by 
1:3) on meshes with N « 10,000 cells; from Becker and Rannacher [30]. 
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Table 3.1: Results obtained by fjL 2 and rjfc ; from Becker and Rannacher [30]. 



10 

nt* 

TOL 

N 

e 

to 

*** 

C s 

N 

M 


4" a 

2836 

6.40-10-' 2 

2.3210- 1 

3.62 

64 

1.47-10 -1 

1.5210 -1 

4” a 

5884 

2.1310 -2 

1.2110" 1 

5.68 

148 

1.0810 -1 

9.80* 10~^ 

4-4 

15736 

7.3610 -3 

4.7610 -2 

6.46 

220 


l&EUiBI 

4 " 5 

23380 

5.5910 _a 

3.1210 -2 

5.58 

592 


2.59-10~ 2 

4"* 





892 

1.19-10“* 


4“ 7 





2368 

5.11* 10 -3 

7.17-10"* 

4~» 

IBHI 




3640 

2.53-10" 3 

3.7210 -3 


Figure 3.4 shows the (scaled) error distribution on meshes obtained by the two 
estimators. The results shown in Table 3.1 indicate that efficient control of the L 2 - 
norm error in the case of heterogeneous coefficients requires the use of ‘weighted’ 
a posteriori error estimates, i.e., the dual weights should be explicitly kept in the 
estimator and evaluated computationally. 

Remark 3.4. ( Curved boundaries) All examples considered so far had been posed 
on polygonal domains which can be exactly matched by the finite element mesh 
domain, i.e., 

flh := U {K G Th} — fh 

This assumption largely simplifies the error analysis and the resulting error esti¬ 
mators. However, in many practical cases at least parts of the boundary of the 
domain are curved and cannot be matched exactly by a polynomial approxima¬ 
tion, e.g., in the cylinder-flow examples presented in the Introduction. Therefore, 
we have to deal with this complication. 




Figure 3.5: Standard situations of cells with curved edges. 
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We consider the two typical situations depicted by Figure 3.5 in which a cell 
K at the boundary has a curved edge Tk C dQ . The computational domain can 
be made to satisfy Qh = ^ , by simply extending ot truncating the domain of defi¬ 
nition of shape functions. Shape functions and transformations are left unchanged. 
This approximation results in a non-conforming finite element scheme, 

(Vw/,, VV’ft) = {f^h) € V h , (3.20) 

in which Vh <jtV, since the elements of Vh will usually not satisfy zero boundary 
conditions. For the error analysis, we assume that the error functional J(-) has an 
L 2 representation, i.e., there is a j e L 2 (Q) , such that J(ip) = (<p, j) . Situations 
in which this is not the case require special considerations. 

Using the solution z E V of the dual problem 

-A z=j in Q, z\ dn = 0, (3.21) 


we obtain the error identity 

J(e) = (j,e) = (e,-A z) = (Ve, Vz) - (e,d n z) dQ . 

For any VVi € Vh , there holds 

(Ve, V^) = (Vu, Vi/> h ) ~ (Vu h ,Vi/>h) 

= (-A u,i/> h ) 4- (d n u,^ h ) dQ - (f.'tph) = {d n u^ h )dn. 


and consequently, 

J(e) = (Ve, V(z-^h)) ~ (d n u,z-^ h )dn + (u h , d n z) dn . 
Now, integrating by parts on each cell K E T h , we obtain 

J(e) = {(/+A%,^-^)k + (d n e,z-il;h)dK} 

Ker h 


(d n U, Z Iph^dQ H" (^/i5 d n Z^QQ . 


This implies the error representation 


J(e) 


E {(Rh, z~^h)K + {rh,z-^ h ) dK } + (uh,d„z) d n, 


KeT h 


with cell and edge residuals defined as above by Rh\k •= f+Aun and 


r h\r • 


\[d n Uh] if r c dK\dQ , 
—d n Uh if T c dQ. 


(3.22) 


We note that the situation considered above is not very practical since the eval¬ 
uation of the discrete equations (3.20) requires the use of numerical integration 
which in general introduces further errors that are not considered here. 
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Remark 3.5. (Nonhomogeneous Dirichlet data) Nonhomogeneous Dirichlet data 

u = g on 

are usually treated by introducing a representative function u G H 1 (Q) satisfying 
u\dn — 9 -> and a corresponding finite element approximation Uh which is the 
interpolation (or the L 2 projection) of g along dQ. The solution is then sought 
as Uh = Uh + u Q h , where u Q h again has zero boundary values. Then, assuming again 
the domain to be polygonal or polyhedral, the error representation (3.8) takes the 
form 


j(e) = {(Rh,z-iph)K + (rh,z—iph)dK) - (g~9h,d n z)dn- (3.23) 

Ke T h 

Remark 3.6. ( Neumann boundary conditions ) The treatment of Neumann bound¬ 
ary conditions 

d n u = g on Tat C dQ , 

does not cause any problems even in the case of a curved boundary, provided 
that the mesh T h is compatible with the decomposition of the boundary dQ = 
T D U Tn • In this case, we simply assume the cells adjacent to the Neumann 
boundary to have possibly a curved edge or face matching the boundary exactly. 
Now, the variational formulation reads 

a(u, ip) = (/, ip) + ( g, ip)dti N € V, 

where V := {v € H l ( f2), W|p D = 0}. For the error of the corresponding Galerkin 
approximation, again Galerkin orthogonality holds, i.e., a(c. v ] h) = 0 for %ph € 
Vh C V. Then, the error representation (3.8) remains valid with the only modifi¬ 
cation that the edge residuals are now defined by 

r ±[d n u h ], if T C dK\dfl, 

Oiir := < 0, if T C r D , 

I g - d n u h , if r C r N . 


3.4 Higher-order finite elements 

We briefly describe the use of the DWR method for higher-order finite elements 

and will see that it can be used in this case without essential changes. Let vj^ C V, 

be finite element spaces of order p+ 1, i.e., they possess the local approximation 

/ \ 

properties of polynomials of degree p. We recall that by setting x/>h = Ih z 
(3.8), we have: 

J(e) = ^2 {( R h,z-I ( h p) z) K + (rh,z-I { h p) z) 9K }, 

K€T h 
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with the cell- and edge-residuals Rh and rh as defined above and some local 

interpolation 1^ z G . In order to extract cell-error indicators for local mesh 
adaptation, we may proceed as before in the low-order case: 

|J(e)| < £ \{R h ,z-I^z) K + {r h ,z-I^z) aK \ 

KeT h 

- X/ {W^WkWz-I^zWk + \\rh\\dK\\z-Ih )z \\dK}- 

K€T h 

This does not reduce the asymptotic efficiency, as will be be demonstrated for the 
special case p = 2. First, we collect the following estimates for the cell residuals 
and weights which are obtained by using standard trace estimates: 

\\Rh \\k = 11/ + A^IIk = ||Ae||/c 
h K 1/2 \\rh\\dK = ^h K 1/2 \\[d n e\\\dK 

< ch~ K 1/2 {h~ K 1/2 ||Ve||* + 4 /2 ||V 2 e||^} 

< c/j^ilVell^ + ||V 2 e||^, 

where the prime in • refers to a cell-wise evaluation of the norm, and K 

1 \ 

is a patch of cells around K . Further there holds the higher-order interpolation 
estimate 

\\z-I { h p) z\\ K + h]( 2 \\z-I^z\\ dK <ch k K \\V k z\\ k , k = 1,2,3. 

Using the above estimates, we obtain on a quasi-uniform mesh, i.e. for hx ~ h , 
that 

|J(e)| < c{*- x ||Ve|| + ||V 2 e|| , }* fc ||V fc z|| < c/i fc+1 ||V 3 u|| ||V fc *||, k = 1,2,3. 

We evaluate this error estimate for the special output functionals 

J(<p) ■■= (V^Ve)HVeir 1 , J{y) := e)||e||- 1 , 

and 

Jipiv)'■= HVV’ir 1 , ^>e/fo( n ). 

which correspond to the energy-norm error, the L 2 -norm error, and the error in 
the H~ l norm || • ||_i , i.e. the norm of the dual space H~ 1 (Q) of Hq(Q) . In 
virtue of the corresponding a priori bounds for the dual solution (for sufficiently 
regular domains Q ), we obtain the L 2 and energy error estimate 

l e ll + /i||Ve|| < ch 3 ||V 3 u||, (3.24) 


as well as the ‘negative-norm’ error estimate 




which are all of optimal order. 
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3.5 Exercises 


Exercise 3.1. Functional-oriented a posteriori error estimates can also be stated 
in terms of energy-norm error bounds. Consider the Poisson model problem. With 
the primal and dual errors e := u—Uh and e* := z—Zh, respectively, there holds 



= (Ve, V*) = (Ve, Ve*) < ||Ve|| ||Ve* 


Then, any a posteriori bound for the energy-norm error supplies also a bound 
for J(e). Specify a situation in which this simple minded approach is inefficient. 
Why is this approach not suited to extract refinement indicators from the error 
estimate? 


Exercise 3.2. Let Q C R 2 be a convex polygonal domain. Develop a residual- 
based a posteriori estimate for the error e := u—Uh with respect to the L°°-norm 
employing a global stability constant. Use the weighted a priori estimate 



< cn |log(e)|, 


* 


for the regularized Green function g ® . 


Exercise 3.3 (Practical exercise). Consider the Poisson problem 

—A u = 1 in f2, u\Qn = 0, 

on the domain Q C (—1, l) 2 , shown in the figure below. Compute the function 
values u(a) and d\u(a) at the point a = (.75, .75). 



a) Use the provided program for computing d\u(a) for a sequence of tolerances 

# 

TOLi = 4“*, i = 1,2,..., and monitor the behavior of the true error and the 
number of cells depending on the achieved accuracy TOL. The code uses the 
duality-based weighted error estimator described in these Lecture Notes and the 
fixed error fraction strategy for mesh refinement. 
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b) Implement your own strategy of mesh refinement based on error indicators, 
smoothness indicators, or a priori information of your choice and try to beat the 
efficiency of the automatic error control process of (a). 



Chapter 4 


Practical Aspects 


In this chapter, we discuss several aspects of the practical use of the DWR method 
described in the previous sections. These are (i) the practical and efficient evalua¬ 
tion of the a posteriori error representations, (ii) the extraction of local refinement 
indicators, and (iii) the design of strategies for economical mesh adaptation. 

The starting point is the a posteriori error representation 

J{e) = {{Rh,z-4>h)K + {rh,z-tl) h ) dK } =: E{u h ). (4.1) 

Ke T h 

as derived above for the Poisson problem. For its evaluation, due to Galerkin or¬ 
thogonality, the subtraction of the arbitrary element ^ £ Vh could be suppressed. 
From this, we extract local ‘refinement indicators’ called t)k of the form 

Vk := | {Rh,z-iph)t< + (rh, z-iph)dK\, 

which are used to steer the mesh adaptation. Collecting these error indicators, we 
obtain an a posteriori error estimator , i.e. an upper bound for the error, 

|J(e)| < 7? := t 1k- ( 4 - 2 ) 

K€T h 

We emphasize that in practice it does not make much sense to estimate further in 
the indicators t]k , since we would inevitably lose sharpness of the error bound. 

For practical use of the error representation (4.1) and the error bound (4.2), 
we have to approximate the terms involving the unknown dual solution z resulting 
in approximate error representations E(uh) and error indicators t)k • This may 
be expensive while the evaluation of the residuals R^ and rh is usually cheap. 
We have to distinguish two related questions: 

• Sharpness of the approximate error representations E(uh) ? 

• Effectivity of the approximate local error indicators t)k for mesh refinement? 
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For measuring the accuracy of the resulting error estimators, we will use the so- 
called (reciprocal) effectivity index defined by 



which represents the degree of overestimation and should desirably be close to one. 

Remark 4.1. Often the physical quantity to be computed can be expressed in 
different forms which coincide on the ‘continuous’ level but differ from each other 
for ‘discrete’ functions and may lead to more or less robust approximations. A 
typical example is the mean normal flux, and, more interesting, the drag and 
lift coefficients computed from solutions of the Navier-Stokes equations. Both can 
be expressed as surface or, after integration by parts, as volume integrals. If the 
desired output functional is not properly defined on the solution space V , such 
as point evaluation in two or more dimensions, it requires regularization, 


J £ (u) = J(u) + o(e), e = TOL. 

In the following, we will mostly suppress the index e indicating this regularization. 

Remark 4.2. In general, the dual problem has to be solved numerically. In the 
case of a nonlinear problem the first step is linearization which will be discussed 
in Chapter 6. Then, this ‘linearized’ dual problem is approximately solved in some 
finite element space Vh C V, 

z h eV h : a'(u h )(tp h ,z h ) = J(tp h ) Vip h €V h , 

where not necessarily Vh = Vh . This will be discussed in more detail below. 

Remark 4.3. If on the basis of a numerical approximation to the dual solution z 
an approximate error representation E(uh) has been generated, one may hope to 
obtain an improved approximation to the target quantity by setting 


J(u h ) := J{u h ) + E(u h ) « J(u). 

This ‘post-processing’ can significantly improve the accuracy in computing J(u ), 
but the resulting error can then no more be estimated on the basis of the available 
information. 


4.1 Evaluation of the error identity and indicators 

Consider the approximation of the model Poisson problem 

-A u = f in 0, u\ dn = 0, (4.3) 

on a polygonal (or polyhedral) domain Q C W* (d = 2 or 3) by piecewise bilinear 
(or trilinear) finite elements (in short ‘Qi elements’). For this prototypical case, 



4.1. Evaluation of the error identity and indicators 


43 


we will compare several common strategies of evaluating the error representation 
(4.1) and the corresponding local error indicators. Notice that the approximation 
of z in E(uh) simply by its Ritz projection Zh £ Vh does not work since, by 
Galerkin orthogonality, it would result in the useless approximation E(uh) = 0. 

Approximation by a higher-order method 

A first possibility is to solve dual problem by using biquadratic finite elements on 

the current mesh yielding an approximation z^ G to z. This yields the 

approximate error representation 

E w (u h ):= ^ {(R h ,z™-I h z™) K + (r h ,z^-I h z^) dK } , 

KZT h 

and the corresponding local error indicators 

Vk* = ( R h,z { h ) -hz { h ) ) K + (r h ,z^ ) -I h z^ ) ) dK . 

For the special situation considered below (see Table 4.1), this approximation 

results in an asymptotically optimal effectivity index, limroL->o = 1* This 
observed behavior is supported by the theoretical analysis presented in the next 
section. However, approximating the dual problem by a higher-order method than 
used for Uh seems not very attractive and can actually be avoided in most cases. A 
compromise would be to compute only the residual with respect to the coarse-grid 

space v!jfo (biquadratics on mesh T 2 h) of the bilinear Ritz projection Zh G Vh 
and to perform a few steps of defect correction using the system matrix of Vh 
with appropriate blockwise preconditioning. 

Approximation by higher-order interpolation 

A simplification is achieved by patch-wise higher-order interpolation of the bilinear 
Ritz projection Zh € Vh of z. Here, we only consider biquadratic interpolation 
in 2-D. On square blocks of four neighboring cells the 9 nodal values of Zh are 

U8€|d to define a biquadratic interpolation 1$ z h • This is then used in the error 
representation instead of z , resulting in the approximate error representation 

E^ 2 \uh) := I^h z h~ z h) k ( r/l ’ ^2h Zh ~ Zh )dK > \ ’ 

KtTh 

and the corresponding local error indicators 

VK ifihi ^2h Z h ~ Zh) ^ {^hi ^2h ~ dK * 

Fbr the special situation considered below (see Table 4.1), this approximation also 

( 2 ) 

results in an almost optimal effectivity index, lim TOL ^ 0 ~ 1. 
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Remark 4.4. In the case of higher-order finite elements , with p > 2 , the evaluation 
of the error identity (4.1) may be done in a similar way as for p = 1 employing a 

patch-wise interpolation ' z h , with p' > p , of the Ritz projection Zh € v fc (p) . 
For example, in the case of biquadratic elements (p = 2), on a 2 x 2-cell patch we 
have 25 nodal values which can be used to construct an interpolation of degree 

p = 4, i.e., in this case we would use z h • 

Approximation by difference quotients 

The error representation (4.1) is estimated by 

\E(uh)\ < ^ Pk“>k, 

KeT h 

with the notation of Proposition 3.1. Applying the usual cell-wise interpolation 
estimates, we have 

u lk ~ \\ z ~ IhA\ 2 K + hK\\z — Ih,z\\dK — c]h)<\\V 2 Z\\ 2 K , 

with an interpolation constant cj ~ 0.1... 1. Now, the second derivatives V 2 z are 
replaced by suitable second-order difference quotients V\zh of the Ritz projection 
Zh of z . This may be even more simplified to the cell-error indicators 

r) { K = cih% 2 p K (u h )\\{d n z h }\\ dK . 

For the corresponding error estimator 

E {3) (u h ) := a ^2 h3 K 2 pK\\[d n Zh\\\dK, 

KeT h 

we usually observe strong over-estimation, i.e. 1^ 1 (see Table 4.1) depending 

on what value we set for the interpolation constant c/ . 

Approximation by local residual problems 

On each cell K , we solve the local Neumann problem (see Bank and Weiser [17]) 

(V%,V^)k = (R/i,V7i)/C + (Th^h)dK V'&h G V/c, 
where Vk — {q G Q 2 (K ), Q-LQi(K)} • Then, in view of the relation 

\{Vv K ^{z-I h z)) K \ < ||V^|| / c||V(^-/^))k|| 

< Cl h K \\Vv h \\ K \\V 2 z\\ K 

w Clh][ 2 \\VVh\\K\\[dnZ}\\dK, 
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the local error indicators may be defined as 

Vk } = Cjh l K 2 \\^V K \\K\\[d n Z h \\\dK- 

In this way, one obtains fairly good bounds for the error in the energy norm. 
However, for local error functionals the results are not better than with the other 
simplified estimators. In such cases the crucial point seems more the accuracy in 
the approximation of the dual solution than in the evaluation of the residuals of 
Uh (see Backes [13]). 


Numerical test 


The effectivity of the first three of these error estimators has been tested for the 2- 
D model problem (4.3) with the solution u(pc) = (1 — x\)(\ — x 2 ) sin(4:ri) sin( 4 :r 2 ) 
for the two output functionals 



Table 4.1 shows the corresponding effectivity indices obtained on sequences of 

locally refined meshes on the basis of the local error indicators 77 ^, i = 1,2,3. 
It turns out that the cheap interpolation-based estimator E^(uh) is almost as 
effective as the more expensive estimator E^(uh)- Therefore, in the following, 
we will almost exclusively use the first one in the presented numerical tests. 


Table 4.1: Effectivity of weighted error indicators for the mean error J\(e) (left) 
and the point-error J 2 (e) (right); from Richter [123]. 
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Remark 4.5. The traditional energy-norm error estimator tje is known to be not 
only reliable , i.e., it provides a safe upper bound for the energy error, but also 
efficient , i.e., it is asymptotically sharp in the sense that 

c 1 ||Ve||<7 7E <c 2 {||Ve|| + ||/-M|}, 
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where fh is the piecewise constant interpolation of / (see, e.g., Verfiirth [132]). 
This means that the estimator is up to a constant asymptotically correct. A cor¬ 
responding result is not possible in general for estimators of locally defined 
error quantities. Already the transition from the error representation to the error 
estimate in terms of (non-negative) cell-wise error indicators is critical, since by 
this localization the asymptotic sharpness of the global error representation may 
get lost. To illustrate this, consider the case J(u) = u(0) and assume that the 
exact as well as the approximate solution are anti-symmetric with respect to the 
aq-axis. Then, e(0) = 0, but usually ^2 Ke j h t)k ^ 0. 


4.2 Mesh adaptation 

Next, we address the practical aspects of successive mesh adaptation. Suppose 
that on the meshes , we have local error indicators rjx extracted from an a 
posteriori error estimate 

|J(e)| < V := J2 VK, N := #{K € T h }. (4.4) 

KeT h 

Using this information the computational mesh may be adapted using various 
different strategies. For quadrilateral meshes, as considered here, the refinement 
and coarsening is facilitated by using hanging nodes . The global conformity of 
the finite element ansatz is preserved since the unknowns at hanging nodes are 
eliminated by interpolation between the neighboring ‘regular’ nodes. 



Figure 4.1: Q\ nodal basis function on a patch of cells with hanging nodes 

We note that there are several alternative strategies for realizing local mesh 
refinement. The occurrence of hanging nodes can be avoided by using special 
‘transition cells’ (triangles or quadrilaterals) which bridge from cells of width h 
to those of width hj2 . The construction of such cells may be complicated in 3-D. 
It may also cause a spreading of the refinement zone which implies extra work 
for de-refining. However, it is basically a question of taste which technique of 
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mesh organization one prefers. Refinement and coarsening of quadrilateral meshes 
involving the use of hanging nodes proceeds as indicated in Figure 4.2. 






Figure 4.2: Refinement and coarsening in quadrilateral meshes. 


At first, we have to check whether on the current mesh the stopping 
criterion 

r)<TOL 

is already satisfied. If this is the case, then uh is accepted as approximation to u 
that represents the target quantity J(u) by J(uh) within the desired tolerance 
TOL . Otherwise, the next refinement cycle is started. To this end, the cells in the 
current mesh are ordered according to 

VKi > • * • > Vk % > * • * > Vk n • 

then, the mesh adaptation may be organized using one of the following strategies. 

Error-balancing strategy 

pelow, we will give an argument which indicates that an optimal mesh is charac¬ 
terized by equilibrated error indicators, such as 

TOL 

VKi ^ } i 1,... A/", 



v( u h) ~ TOL. However, the so-called error-balancing strategy based 
H^this idea is ‘implicit’ as it involves the current number of mesh cells N which 
obtained only at the end of the adaptation cycle. We check, starting from i = 
j = 0, whether 

TOL 

Vk 




< 


N + 3j 
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If this is not satisfied, then the cell is refined, the counters j and i are 
increased by one, and one proceeds to the next smaller 77 /^ . But, if the condition is 
satisfied, then the new mesh J™ ew is reached. This strategy is potentially optimal 
but involves many expensive checking operations and is therefore impracticable. 


Fixed-error-reduction or fixed-rate strategy 

For fractions X , Y, with 1 — X > Y, determine indices N*,N* G {1 ,...,7V}, 
such that 

N* N 




Xn, 


E 11 k 


Yn. 


i= 1 


i=N, 


Then, the cells K \,..., Kjy* are refined and the cells ,..., Kn are coarsened. 
Common choices are X = 0.2 and Y = 0.1. Alternatively, in the fixed-rate 
strategy , one refines X-N and coarsens Y-N cells with largest and smallest error 
indicators, respectively. For appropriate choices of X, Y this accomplishes to keep 
the number of cells almost constant in the course of the mesh adaptation process. 


Mesh-optimization strategy 

The information contained in the error representation (4.1) may be used directly to 
construct an ‘optimal’ mesh on which the error tolerance ry « TOL is achieved, i.e. 
skipping the intermediate ‘one-level’ refinement steps. Here, a mesh T h = {K} is 
characterized by a continuous mesh-size function h(x ), where h\x « hx . Further, 
it is assumed that the error estimator is related to a continuous limit of the form 


E 

KeJ h 


Vk 


/ h(x) 2 $(x) dx =: 77 (h), 
Jn 


with $ := (<FiH-^> 2)^ > 3 a mesh-independent weighting function. We note that the 
latter requirement rules out cases in which the functional J(-) and hence also the 
dual solution implicitly depend on the mesh size, as for example in the estimation 
of the error in the energy or the L 2 norm (cf. Section 3.2). The components of 4> 
are defined by the limiting processes of residuals and weights, for TOL —► 0: 


max \ f+Auh\ 


max $1 

xeK 


\h K l max|[ d n u h ] 9K 


max $ 0 , 

x£K 


h J 2 max 

K xeK 


z-hA 


max $3 

xeK 


(4.5) 

(4.6) 


Here, the mesh-size power h(x) 2 is related to the ‘order’ of the considered finite 
element ansatz. The justification of these assumptions will be discussed in the next 
chapter. The function <£(#) may have strong (regularized) singularities which will 
possibly require the mesh-size h(x) to reduce down towards zero even for tolerance 
TOL > 0. Further, we introduce the mesh complexity formula 


N 


E h 


d h d 
K a K 


KeT h 


[ h (x) 

Jn 


-d 


dx =: N(h), 
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which denotes the number of degrees of freedom for a mesh-size function h(x). 
With this notation, we obtain the following result. 

Proposition 4.6. The mesh-optimization problem 


r)(h) —> min, N(h) < Af, 


max 


is solved by 


bopt (%) 


W \l/<* 


N, 


*<*) 


-1/(2+d) 




max 


provided that 


W := f $(x) d/(2+d) dx 

Jn 


< oo. 


For an 1 optimal 1 mesh, there hold 


VP( 2 +d)/d j^(2+d)/2 

TOL = —-T7TT-;— and N 


N 2 / d 


TOL d / 2 * 


(4.7) 


(4.8) 


(4.9) 


Proof. Following the classical Lagrange approach, we introduce the Lagrangian 

L(h, A) = r)(h) + \{N(h) - JV max } 


with Lagrangian multiplier A G M. Then, the optimal mesh-size function h opt is 
characterized as stationary point of the first-order optimality condition 


d_ 

dt 


L(h-\-tip, A)|* = o — 0? 


d_ 

dt 


L(h, A+t/x)|i = o — 0, 


for all admissible variations and p. This means that 


2 h(x)$(x) — d\h(x) 


-d-1 


0, 


[ K*) 

Jn 


-d 


dx — TV, 


max 


o, 


and, consequently, 


h(x) 


2 


dX 


$(*)) 


— l/(2+d) 


2 \ d / (2+d) 


* 


d\ 


f $(x) d/(2+d) 

Jn 


dx = AT, 




Rrom this, we deduce the desired relations 


A = - h(x) 2+d *(x), 


b*opt (^r) 


W \ i/rf 


N, 


$(x) 


—l/(2+d) 


max 


Rrom the formula for /i opt , we conclude that on the optimal mesh, there holds 


TOL 


W \ 2/d 

~N ) 


f $(x) 2 ^ 2+d ^$(x) dx 

Jn 


W<2+d)/d 

m/d 


1 



which proves (4.9). 
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Remark 4.7. We note that in two dimensions the optimal mesh complexity is 

TOL = ©(AT 1 ) or N = 0{T0L~ l ), 

for linear or bilinear finite elements, provided that sup TOL ^ 0 W < oo is satisfied. 
This is the case even for rather ‘irregular’ functionals J(-). For example, the eval¬ 
uation of J(u) = diu(a) leads to $>(x) « {\x — a\ 2 + TOL 2 )~ 3 and, consequently, 

sup W « / \x — a\~ 3 / 2 dx < oc. 

TOL—>0 A* 

In the case that sup TOL W = oo , as for example for J(u) = d 2 u(a) , the optimal 
mesh complexity becomes 

TOL = 0(N a ~ 1 ) or N = 0(T0L~ a ~ l ), 


with some a > 0. In particular, for the latter functional, we easily find TOL = 

0{N-'\\og{N)\). 

Remark 4.8. The result of Proposition 4.6 implies that the optimal mesh-size 
distribution is characterized by the equilibration property 

h{x) 2+d ${x) = (^^) (2+d)ld = const. (4.10) 

This justifies the strategy of equilibrating the local error indicators t)k , 

T)(h) = ( h(x) 2 $(x) dx « ^2 h2 K d $K= Y! 

Jn KeT h Ke Th 

as used in the error-balancing strategy. 

Let the weight function $(#) be bounded. Then, once a balanced mesh 
satisfying (4.10) is reached, a maximum increase of accuracy is achieved by uniform 
mesh refinement. If & is singular, then the ‘optimal’ mesh size tends to zero at 
these singularities even for TOL > 0. 

Remark 4.9. Alternatively to (4.7), we can also consider the mesh-optimization 
problem 

N(h) —> min, rj(h) < TOL, 

which has the solution 
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Remark 4.10. Although the mesh-optimization strategy seems very attractive, its 
realization involves several problems: 

• The derivation of the formula for an optimal mesh-size distribution is based 
on the assumption that on the considered meshes the cell-residuals behave in 
an optimal way under refinement, i.e., pk ~ hx . This is hard to prove and 
may not be true in general; see the discussion of this question in Chapter 5. 

• The numerical approximation of the weighting function $(x) should provide 
more information than can be cheaply obtained using only information on 
the current mesh. 

• The explicit formulas for h opt (x) have to be used with care in designing a 
mesh as their derivation implicitly assumes that they actually correspond to 
scalar mesh-size functions of isotropic meshes, a condition, however, which 
is not incorporated into the formulation of the mesh-optimization problems. 


A detailed study of how to utilize mesh optimization for the model problem has 
been made by Richter [123]. 

Remark 4.11. (The role of regularization) At the beginning of this chapter, we 
had mentioned that in the case of a ‘singular’ functional like, for example, that 
for evaluation of the derivative point value, regularization is necessary, 





d\udx, 


To demonstrate the effect caused by regularizing with e = TOL, as proposed 
above, compared to e = h m i n , we show the results of a numerical test for the 
model situation considered in Example 3.4, see Table 4.2. We see that in this case 
regularization drastically reduces the number L of refinement steps for reaching 
the same accuracy level which means less computational work. 


Table 4.2: Results for computing <9iu(0) in Example 3.4 using regularization with 
e = TOL (right) and £ = hmin (left); from Becker and Rannacher [30]. 



£ — hmin 

e = TOL 

TOL 

N 

L 

die(0) 

N 

L 

<9ie(0) 

1 

40 

4 

2.57-10 -0 

40 

4 

2.57-10~° 

4" 1 

124 

6 

7.3810 -1 

64 

4 

1.47-10-° 

4~ 2 

448 

11 

1.3510 -3 

148 

6 

7.5110 -1 

4 6 

1780 

15 

1.19-10 -3 

940 

9 

4.1010 -1 

4 -4 

6328 

19 

2 .20-10 -3 

4912 

12 

4.1410 -3 

4“ 5 

25984 

24 

4.28-10 -4 

20980 

15 

2.27-10 -4 

4" G 

95260 

28 

1.39-10 -4 

86740 

17 

5.82-10 -6 







































52 


Chapter 4. Practical Aspects 


4.3 Use of error estimators for post-processing 

The duality approach described so far for estimating the discretization error can 
also be used for increasing the accuracy in the approximation of the target quan¬ 
tity. Consider the model situation as before, i.e. the Poisson problem in 2-D written 
in variational form as 


a(u, V>) = (/, Ip) ^ € V, (4.11) 

and discretized by 

a(U' h ^ h ) = {f,'ip h ) ip h e V h . (4.12) 

The output functional is J(-). Let z € V be the corresponding dual solution and 
Zh G Vh its Ritz projection on the ‘primal’ mesh . By definition, there holds 

J(u) = a{u,z) = (/,*), (4.13) 

1. e., if z were known, the target quantity J(u) is determined solely by the data /. 
This observation will be used below for post-processing the approximation J(uh) • 
For this, we recall the identity 

J{u) = J(u h ) + a(e,z) = J(u h ) + p{u h )(z-z h ), 

where 

p(u h )(z-z h ) = if,z-z h ) - a{u h ,z-z h ). 

Above, we have discussed ways of constructing approximations to the dual solution 

( 2 ) 

2 , e.g. the patch-wise biquadratic interpolation Zh := I^h z h of Zh on the mesh 
Tfc . This led us to the approximate error representation 

J(e) « p(uh)(zh-z h ) = p{u h )(z h ). 

Rewriting this relation as 

J(u) SS J^Uh) := J(u h ) + (/, Zh) - a(u h ,z h ), 

we obtain a presumably better approximation to J(u) than is ). The error 
can be written as 

J(u) - Ji(u h ) = J(e) - p{u h )(z h ) = a(e,z) - a{u h ,z h ) = a(e,z-z h ), 
which implies 

\J{u) - AH)! < l(Ve(( !!V( Z -4)||. (4.14) 

Since it is not clear whether Zh is a reasonably better approximation to 2 than Zh , 
this estimate is of only questionable value. Further, the two energy-norm errors 
correspond both to the ‘primal’ mesh T h and can therefore not be minimized 
independently. It would be desirable to have the possibility of using independent 
meshes for Uh and Tj£ for Zh in constructing an approximation of J(u ). 
This is accomplished in the following proposition (see Richter [123]). 
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Proposition 4.12. Let T h and T£ be two independent meshes and Vh and Vj* 
corresponding finite element spaces in which the Ritz projections Uh and z £ of u 
and z are computed. Farther, denote by z £ the patch-wise biquadratic interpola¬ 
tion of z^ on the dual mesh T£ . Then , for the post-processed approximation 

hM := J(u h ) + (f,z* h )-a(u h ,z* h ), 


there holds the estimate 

\J(u) - J 2 (u h )I < ||Ve|| \\V(z-z* h \\. (4.15) 


Proof. We have 

J(u) - J 2 (u h ) = J(u)~ J(u h ) - p(u h ){zD 

= (/, z) - a(u h , z) - (/, z* h ) + a(u h , z* h ) 
= {f,z~Zh) ~ a(u h ,z-zl) 

= a(e, z — z£). 


This implies the assertion. □ 

We emphasize that in the estimate (4.15) the two energy-norm terms can be 
minimized independently by optimizing the primal and dual meshes and T£ , 
respectively. One may hope to obtain an even better approximation by 

M u h) := Muh) = J(uh) + (f,Zh) ~o-(uh,zl), 

where Uh is the patch-wise biquadratic interpolation of Uh . Below, we will test 
all three post-processed approximations to J(u) for the model problem on Q = 
(—1,1) 2 with the prescribed solution u(x) = (1— x\)(\—x^) exp(l— x^ 4 ) and the 
output functional 

J(u) = L u{x\, 0.5) dx. 

In this example primal and dual solution have irregularities at different locations 
such that it is to expected that maximal efficiency is achieved using different 
meshes for Uh and Zh as shown in Figure 4.3. 

Figure 4.4 shows the mesh efficiencies of J{uh) and the three post-processed 
approximations J\(uh ), J 2 (uh) , and J${uh ). The mesh refinements are driven 
by energy-norm error indicators as derived in Section 3.2 separately on the primal 
and dual meshes and T£, respectively. We see that J\(uh) indeed does not 
bring significant advantages over the original approximation J(uh ). The two other 
approximations J 2 {uh) and Jz(uh) which use different meshes for u and z are 
clearly superior and show a mesh complexity like TOL « N~ 2 (N = N pr i ma j + 
^Vduai)- This is to be compared to the optimal complexity TOL « N~ l obtained 
above in Proposition 4.12. 
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Figure 4.3: Primal (left) and dual (right) solution of the model problem; from 
Richter [123], 



Nprimal + Ndual 


Figure 4.4: Mesh efficiencies of post-processing: J(u/ l ) (solid line, symbol +), 
J\{uh) (broken line, symbol *)> (broken line, symbol *), Jz(uh) (dotted 

line, symbol D); from Richter [123]. 
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4.4 Towards anisotropic mesh adaptation 

Sometimes isotropic mesh refinement as discussed so far is not efficient for properly 
resolving direction-dependent features of the solution. For example, in singularly 
perturbed problems of the form 

-eAu + bu = / in u\ dn = 0, 

with small coefficient e , boundary layers may occur in which the solution has a 
large derivative in normal direction to the boundary while it varies only slowly in 
tangential direction. A similar phenomenon occurs in 3D along reentrant edges of 
the domain (edge singularities). In such a situation it is appropriate to use meshes 
which are anisotropically refined in the sense that the cells along the boundary 
are much thinner in normal than in tangential direction, see Figure 4.5. 



Figure 4.5: Locally anisotropic tens or-product meshes. 


The questions are now whether the weighted error estimator contains information 
about anisotropy in the exact solution, and whether we can extract local indicators 
which tell us how to adapt the mesh according to (i) orientation of cells, and (ii) 
optimal stretching of cells. This aspect of automatic mesh adaptation is a rather 
difficult one and still the subject of current research. 

To fix ideas, consider a given function u € C 2 (B\) on the ball B\ = {x € 
R d , |x| < 1} and determine the direction of smallest line-wise L 2 error of linear 
interpolation, i.e., find the unit vector e such that 







min. 
r 


Here, T = {x G B\,x = se,0 < s < 1}, and H(u) := V 2 u is the (symmetric) 
Hessian matrix of u . Hence, the interpolation error becomes minimal for e = e m i n 
being the eigenvector corresponding to the eigenvalue of H(u) with minimal mod¬ 
ulus. We see that all the information about the best orientation of cells for mini¬ 
mizing the interpolation error is available by the Hessian matrix. The eigenvalue 
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quotient o = |A max /A m i n | determines the cell aspect ratio and the corresponding 
(orthogonal) eigenvectors e m i n and e max the orientation of the cell. 



K 

hi 


Figure 4.6: Cell of a Cartesian mesh . 


Quadrilateral meshes are not very suited for dynamical cell reorientation. For this 
purpose, triangular meshes are more flexible. Therefore, we will consider adap¬ 
tive cell stretching alone and, for simplicity, will concentrate on the construction 
of ‘optimal’ Cartesian meshes consisting of cells K with area \K\ = / 11 / 12 , as 
shown in Figure 4.6. Again ‘hanging’ nodes are allowed. The maximum gk •= 
max{/ii//i 2 , h 2 /h\} is the cell aspect ratio and ah := m&XKej h &k the maximum 
mesh aspect ratio. 

We recall the a posteriori error representation 


J(e) = {(f + & u h, z ~Ih z )K + l([d n u h ],z-I h z) aK \ d( i}. (4.16) 

K€Th 


Practical experience and theoretical analysis show that on isotropic meshes the 
edge-residual terms ([d n Uh], z— Ih,z)dK can be made to dominate the cell-residual 
terms (/ 4- z — Ihz)K , see Exercise 5.2. Let us assume that the same is true 
also for anisotropic meshes (for some theoretical support see Kunert and Verfiirth 
[99]). Therefore, we only consider the edge-residual terms and as before suppose 
the approximation 

h~ l [diUh] ~ d^u {% = 1 ,..., d). 

Then, assuming again the second-order derivatives of u as constant, we obtain 


( 7 2 J-hZfdK 



Minimizing this with respect to h\ yields the necessary condition 


2h\\d2u\ \dfz\ — 2\K\ 2 h 1 3 |^iu| | d%z 





\d\u\ | d\z\ 
\d%u\ \dfz\ ’ 


and, consequently, 


hi ~ ( \dfu\ [ d\z 

h 2 ~ \|^2 U I \d\ z 




1/2 


(4.17) 
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This result is counterintuitive as it does not indicate the optimal cell stretching. 
To see this, consider the case that u is linear in xi-direction, i.e., d\u = 0, and 
that z is isotropic. Then, formula (4.17) would suggest to refine the cell in x\ 
direction which is evidently the wrong decision. It seems that considering only the 
edge terms in (4.16) leads to contradictory results, and we rather have to take into 
account the whole combination of cell and edge residuals. 

Remark 4.13. The development of a rigorous criterion for anisotropic mesh re¬ 
finement on the basis of ‘goal-oriented’ error representations such as (4.16) must 
be left as an open problem. Here the difficulty is caused by the local interplay 
of the primal and dual solutions which may have significanty different regularity 
properties. In the special case of error estimation with respect to the energy-norm 
the dual solution coincides with the primal error such that both quantities behave 
very much the same way. Then, anisotropic mesh refinement can be surely guided 
by information from the jump residuals of Uh alone (see, e.g., Siebert [126]). 

In view of the above discussion, we now follow a more heuristic approach and 
base the anisotropic cell adaptation on an estimate for the interpolation error. We 
recall the anisotropic interpolation error (see, e.g., Becker [18]) 

\\V (u -1 h u)\\ K < c(fc?||ftVu|ft + h2||&Vu|&) 1/2 . (4.18) 

Hence, assuming the second-order derivatives as constant on K , we have 

||V(«-7 fc tt)||jc < c\K\V 2 (hlldxVuf + \K\ 2 hi 2 \d 2 Vu\ 2 ) 1/2 . 

Minimizing this with respect to h\ results in the necessary condition 

2h 1 \d 1 Vu\ 2 -2\K\ 2 hi 3 \d 2 Vu\ 2 = 0 =► h\ = 

and, consequently, 

hi \d 2 Vu\ 
h 2 ^ |3iVu|' 

In view of this result, we now consider the heuristic error indicator 

m ■= ||V(u-4u)||k-||V(z-4z)||k-, 


which is minimized for 


hi _ \d 2 Vu\\ \d 2 Vz 
h 2 ~ |3iVu||3iVz| 



This relation simultaneously reflects possible anisotropies in the primal and dual 
solution. The performance of anisotropic mesh adaptation based on the relations 
(4.17) and (4.19) will be compared in some numerical tests below. 
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Remark 4.14. We note that for the goal-oriented adaptation of purely tensor- 
product meshes an analogue of the mesh optimization strategy described in Section 
4.2 can be developed, see Richter [123]. This approach yields meshes of optimal 
complexity in the case that the anisotropies of the primal and the dual solution 
are aligned with the coordinate directions. An extreme case occurs in Example 
3.5, on a square domain, where the functional 

J(u) = / d n uds 

Jdn 

results in a dual solution which is concentrated along the boundary dCl . Then, an 
optimal anisotropic refinement yields a mesh complexity of the order 

TOL = OK 1 ), 

for constant number N of cells. This is much better than the complexity TOL = 
0(N~ l ) , which can be achieved with isotropic meshes. 


Numerical test 

As test case, we consider the Poisson problem on the square domain Cl = (—1, l) 2 , 

-An = f in fi, u\ dQ = 0, (4.20) 

for the exact solution 


u(x) = {l-x\) 2 (\-X 2 ) 2 (kx\+t)A) x , 

where the parameter k = 1,4,16,64,..., determines the strength of the aniso¬ 
tropy. The right hand side is determined as / := —A u . This solution is shown for 
k = 4 and k = 64 in Figure 4.7. The quantity to be computed is the mean value 

J(u) := / udx. 

Jn 

In this case the anisotropy is only in the primal solution while the dual solution 
satisfies —A z = 1 and is isotropic. 

The computation starts from a coarse uniform tensor-product mesh which is 
then successively adapted on the basis of the relations (4.17) (‘Strategy F) and 
(4.19) (‘Strategy II’). The efficiency of the resulting meshes is depicted in Figure 
4.8. Apparently, the meshes produced by Strategy II on the basis of the heuristic 
relation (4.19) are more efficient as the anisotropy increases than those obtained 
by Strategy I on the basis of the ‘rigorous’ relation (4.17). 



4.4. Towards anisotropic mesh adaptation 


59 



Figure 4.7: Anisotropic solutions for k = 4 (left) and k = 64 (right); from Richter 
[123]. 





Figure 4.8: Mesh efficiencies of anisotropic refinement by Strategy I and Strategy 
Hy for k = 1 (left) and k = 64 (right); from Richter [123]. 
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4.5 Exercises 

Exercise 4.1. The fixed rate strategy in a posteriori mesh adaptation refines and 
coarsens certain fractions X and Y of those cells with largest and smallest indi¬ 
cator values t]k , respectively, 

VKi > * • * > VKi > * * • > Vk n • 

Design refining/coarsening strategies by specifying X , Y, such that 

a) the number N of cells approximately doubles in each refinement cycle; 

b) the number N of cells is approximately kept constant during the refinement 
process. 

Exercise 4.2. The well-known Bramble-Hilbert theory guarantees the cell-wise 
error estimate 

\\V(u-I h u)\\ K < Cl h K \\V 2 u\\ K 

for the piecewise bilinear finite element nodal interpolation I^u on regular meshes, 
with a constant c/ independent of K . Sketch in 2D the argument that this esti¬ 
mate has an analogue on cells with one or two hanging nodes , 

\\V(u — Ihu)\\K < cih K \\^ 2 u\\ k , 

where K is the mother cell of K on the preceding refinement level. Recall that 
at a hanging node the function value is set to be the average of the values at the 
two neighboring ‘regular’ nodes. 

Exercise 4.3 (Practical exercise). Hanging nodes are commonly used to ease mesh 
refinement and coarsening in triangular or quadrilateral meshes. However, the pres¬ 
ence of hanging nodes (and also transition cells) usually destroys the uniformity 
pattern of the mesh which may drastically reduce the accuracy of approximation. 
To demonstrate this, consider the model Poisson problem 

-A u = f in D, u\ d n = 0, 

on Q = (—1, l) 2 , with right-hand side f(x) = \n: 2 cos(^Xi7r) cos(^X27t) and 
corresponding exact solution u(x) = cos(|xi7r) cos(^X27t) . Compute the point- 
value diit(0.5,0.5) on a sequence of 

a) uniform meshes with mesh-sizes /i = 2“ 1 ,i=l,2,...; 

b) locally refined meshes using the a posteriori error estimator and the corre¬ 
sponding refinement strategy of Exercise 2.3. 

Compare the achieved accuracy in terms of the number of cells N . 



Chapter 5 

The Limits of Theoretical 
Analysis 


In this chapter, we want to discuss some questions concerning the theoretical jus¬ 
tification of the DWR method for goal-oriented mesh adaptivity presented so far. 
We will see that this task is rather demanding and poses several new questions for 
the theoretical analysis of the finite element method. In fact, relying on the avail¬ 
able results from the literature, we do not reach very far, yet. Since several not very 
practical assumptions will be used, we dispense with stating formal propositions. 

To illustrate the problem, let us consider the special case of the evaluation 
of the derivative of a smooth solution u € C 2 (fi) at some point a € Cl C R 2 . For 
this, we use the regularized output functional 




diu(a) dx = d\u(a) + 0(s 2 ), 



TOL. 


The corresponding dual solution behaves like a regularized ‘derivative Green func¬ 
tion’ of the Laplacian: 


|V fc z(x)| wrtx)-* -1 



2 + £ 2)-(k+l)/2 


Then, for bilinear elements, the a posteriori error estimate takes the form 


|J £ (e)| « Tf u : 


T, 


Pk , 


Ker h 


rx := maxr(x) 

x€K v y 


Now, assume that the local residuals are related to the local mesh size like 


Pk ~ hx , 


( 5 . 1 ) 
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uniformly for every K G and h > 0. This condition may be checked a poste¬ 
riori in the course of the mesh adaptation process. Then, we obtain 

Y, ( 5 - 2 ) 

KeT h K 


In Section 4.2, we have seen that the optimal mesh for prescribed accuracy TOL 
is characterized by the equilibration property 




E 

Ke T h 


TOL 

N 


TOL. 


From this, we derive 


and consequently, 


h K 


3/2 (TOL\ 1 / 2 




K 


V 


N 




N 


h K h 


-2 

K 


KeT h 


N \l/2 

tol) 


E 

KeT h 


2 
K 


hlsV 


- 3/2 

K 



This implies that 


N oc TOL -1 , (5.3) 

which is better than the N oc TOL ~ 2 that could be achieved on uniformly refined 
meshes on the basis of the general a priori convergence estimate 

\die{a)\ = 0{h). 


This predicted asymptotic behavior is well confirmed by the results shown in 
Table 5.1. We emphasize that in this example strong mesh refinement occurs, 
although the solution is smooth. In fact, this phenomenon should rather be inter¬ 
preted as ‘mesh coarsening’ away from the point of evaluation. Further, observing 
that r m in ~ TOL and r max ~ 1, for the optimized mesh, there holds 


h 


min 


TOL 5 / 4 , /i max » TOL 1 / 2 , 


and, consequently, h mm 


w J , 5 /2 

~ '‘'max > 


which means that 


5 log (TOL) 
4 log(2) 


refinement cycles are needed to reach the ‘optimal’ mesh. 
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Table 5.1: Computing d\u(0) using the estimator r (L refinement levels); from 
Becker and Rannacher [30]. 


TOL 

N 

L 

Me) 


4-3 

940 

9 

4.1010 -1 

1.42-10 -2 

4 - 4 

4912 

12 

4.1410 -3 

3.5010 -3 

4-5 

20980 

15 

2.2710' 4 

9.2510- 4 

4-5 

86740 

17 

5.82-10 -8 

2.3810- 4 


Remark 5.1. Relation (5.1) is decisive for the optimality of a refined mesh. To see 
this, suppose that only 


PK 


h l r 


holds, for some small e > 0. Then, the above calculation would result in 


^2 _ ^ 6 /( 4 - e )/ < TOZ /^ 2 /( 4 £ ) 


K 


K 




N 




and consequently, 


N 


Y h K h 


-2 

K 


(JL-\ 2/(4 £) V h 2 r _6/(4_e) oc (— _1 

\TOl) 2^, k k \tol) 


KeT h 


N \ 2/(4—e) 


K€ T h 


This would give us the asymptotic complexity 

N oc TOL -W 2 , 

which growths faster than TOL -1 . 

Remark 5.2. An alternative, more explicit strategy for mesh adaptation may be 
based on the balancing condition 


hKPK TOL 

-o— (X 


K 


|fi| 


TOT 

| J e {e)\ oc Y, h 2 K ^rrr- » TOL. 

K<iT h 


ini 


Here, the complexity analysis gives us (assuming again that pk ~ hx ) 


/TOL\!/2 3/2 


hx a 


K » 


and, consequently, 


N 


Y h K h 


-2 M 

K OC - 


K€f h 


TOL 


V h 2 r" 3 oc M 
Ke T h iUL 


2 » 


which shows that this approach is not efficient for very singular error functionals 
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5.1 Convergence of residuals 

The above example illustrates that the assumed asymptotic relation (5.1) for the 
residuals pk is crucial for deriving optimal complexity results. For analyzing 
this question in the case of d-linear finite elements, it suffices to consider the 
edge-residual part ||[d n ^]||a/f» since on Cartesian meshes the cell-residual term 
automatically satisfies 

\\f + Allh\\K = ll/ll K < c{f)hK, 

for bounded /. Now, notice that 

h^ ,2 \\[d n u h }\\ aK =: \D 2 h u h \ K \ 

can be viewed as a mean value of a second-order difference quotient of Uh on the 
cell-patch K containing K and its neighbors. Hence, in order to establish the 
bound (5.1), we have to seek for an estimate of the form 

sup max \D\uh\ < c(u), (5.4) 

h >o 

where the constant c(u) depends on bounds for higher-order derivatives of the 
solution u. For proving (5.4), one may use the local inverse relation, 

||V<?||/r < c r h^\\q\\ K , q e P r (K), 
and the natural nodal interpolation I^u € 14 , 

\Dlu h \ K \ < h^\\D\u h \\ K 

< h^\\D 2 (u h -I h u)\\ K + h^\\D 2 h I h u\\ K 

< ch~ K 2 \\Ve\\ K +ch- K 2 \\V{u-I h u)\\ K + h~ K l \\D 2 h I h u\\ K 

< ch^WVeWK+ch^W^uWn 

< c||V 2 u||oo, 

where again K denotes a cell-patch neighborhood of K . Unfortunately, this argu¬ 
ment only works on quasi-uniform meshes, for which h ma , x /h m [ n < c is assumed, 
since the local error estimate 

l|Ve|| oo\K < hKC(u) 

does not hold in this strong form on meshes with h max /h m [ n —> oo . What one can 
prove is the weaker version 

UVelU* < c(u){ max h K '+h 2 }, 

K’eS(K) 

A 

where h := /i max , and S(K) is some (9(l)-neighborhood of the cell K . Alter¬ 
natives may be found by adopting weighted-norm techniques from the ‘classical’ 
pointwise error analysis of finite element approximation. However, most of these 
results also require the mesh to be quasi-uniform. Hence, the problem must be left 
open. 
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5.2 Approximation of weights 

Next, we want to analyze the effect of approximating the dual solution 2 in the 
weights ujk > on the accuracy of the error estimate. For simplicity, we again restrict 
our attention to the two-dimensional case. In the following, we will examine two 
different methods of approximating z: 

i) Approximation by a higher-order method. First, we consider the approximation 

of the dual solution z by its Ritz projection into the space V^ of biquadratic 

finite elements, defined by 

(V^,v4 2) ) = JWh) ¥> h ev k (2) . 

The resulting approximate error representation then reads 

E{u h ):= E {(• r /»-4 2) - / /»4 2) )^ + (r^z^-hz^^dK}. 

JfeT h 

Its difference to the exact error representation E(uh) can be written in the form 


E(u h ) - E(u h ) = p(u h )(e*-I h e*) = (Ve, V(e*-/ fc e*)), 

with the abbreviation e* := z—z^ . For estimating this error, we recall the well- 
known a priori error estimate for biquadratic finite elements: 


(E ii yfce '* 

K€ T h 



* = 0 , 1 , 2 . 


Using this, we can then estimate as follows: 


E(uh) - E(uh )| < || Ve|| || V(e*-I h e 


<ch\\V 2 u\\( E h K l|V 2 e* 


i) 


1/2 


KeT h 


<ch 3 \\V 2 u\\ ||V 3 z 



This estimate is unsatisfactory, as it requires the primal as well as the dual solution 
to be smooth which rules out most interesting applications. However, this objection 
can be somewhat weakened by the following modifications of the argument: 


E(u h ) - E(uh) 


= |(Ve,Ve*)| = 

<ch 2 ( E h K 

JfeT h 


l(V(u 

IV 2 ull 


hu), Ve*) | 


*) 


1/2 


V 3 z 
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or, setting 



u—u 


( 2 ) 


the ‘biquadratic’ Ritz projection error, 


E(u h ) - E(u h )| = |(Ve, Ve*)| = |(Ve, V(z-I h z)) 


< ch 2 \W 3 u 


( £ k K 

Ke Th 




Here only smoothness of one of the solutions u and z is required and singularities 
in the other one can be compensated by proper mesh refinement. Though the 
estimates (5.6) and (5.7) are better than (5.5), they still do not cover point-error 
evaluation, i.e., J(u) = u(x o), since in this case 


£ 4liv 2 z 


2 

K 


£ h K r K 


0(1), 


KeT h 


K€T h 


and generally h 2 is not smaller than TOL . To get a better result, we may use a 
L°°-L l -duality argument as follows: 


I E(u h ) - E(u h )| = |(Ve, Ve*)| = |(V(u-/f«), Ve*))| 


< cmax{/i| ( -||V 3 u|| oo; /c} f |V 

K Jn 

< ch\\og(h m in)\ max{h 2 K \\V 3 u 

fC 


e* dx 


(5.8) 


oo 


\k}- 


The L 1 -error estimate for the regularized Green function, 



used in deriving (5.8) has been proven in Frehse and Rannacher [60]) only for quasi¬ 
uniform meshes; however, its extension to locally refined meshes is a technical 
exercise. The estimates (5.5)-(5.8) are useful provided that h 3 TOL on the 
current mesh. According to the discussion at the beginning, this is satisfied even 
in the case of derivative point value evaluation with 


TOL 



2 

max' 


However, since computing the dual solution by higher-order elements is too ex¬ 
pensive in most practical situations, we will not pursue this discussion further. 


ii) Approximation by higher-order interpolation. Next, we consider the theoretical 
justification of the approximate error estimator 

E(u h ):= £ {(Rh,l 2 h z h~Zh)K + ( r hJ 2 h )z h-Zh)dK}, 

Kdt h 
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where I^ z h is the patch-wise biquadratic interpolation of the bilinear Ritz pro¬ 
jection Zh as defined above. This raises the question: Why should 4h z h te a 

better approximation to z than Zh ? In fact, the construction of z h is based on 
nodal point information of Zh , and the point error (z — Zh)(a) behaves generally 
not better then (D(h 2 ), even on uniform meshes. Hence, it seems unlikely that 



K < 


Z-Z h 



However, this may not be the right point of view. We could also seek to prove the 
somewhat weaker relation 

\p(u b )(z-l!$z h )\ < \p(u h )(z-z h )\, 


which has the flavor of a global ‘super-approximation’ property. Therefore, it will 
probably depend on some uniformity property of the mesh. In order to pursue this 
thought further, we rewrite the error identity J(e) = p(uh)(z — Zh) in the form 

J(e) = p(u h )(z-l!$z) + p(u h )(l!$z-l!$z h ) + p{u h )(4h z h -z h ), (5.9) 


where the last term on the right is just our approximate error estimator. The first 
and second term will be estimated separately. To this end, we have to assume that 
the meshes have been optimized such that 


Pk < c/ix, K £ T h- 



( 2 ) 

Using the local approximation properties of the interpolation z , 


(I \ z ~4 2 h z 


k + hhieWz—1$ z||Ik) ^ ^ c / 2 ^!dl^ 73z lls(K’)) 


where S(K) denotes the cell-patch on which I^z is defined, we obtain for the 
first term in (5.9): 


\p(u h )(z-1!$Z)\ < c( h Kp 2 K ) 

K€T h 




Consequently, using (5.10) we arrive at the estimate 

\p(u h )(z-I$z)\ <c(u,z)h 3 . (5.11) 

The second term in (5.9) is the ‘hard’ one which requires more work, as it relates 
properties of the non-local Ritz projection and the local interpolation. Its esti¬ 
mation strongly relies on uniformity properties of the mesh . The idea is that 
the scaled error h~ 2 e is a ‘smooth’ function, such that it can be approximated 
in Vh . To make this concept clear, suppose that the mesh Tfc is uniform with 
mesh-width h . Then, it is known that in the nodal points, the error z—Zh allows 
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an asymptotic expansion in powers of h which can be expressed in the form (see 
Blum and Rannacher [34]) 

IhZ-z h = I h (z-z h ) = h 2 I h w + h 3 r k , 

with some /i-independent function w G Hq(Q) and a remainder satisfying Ht/^I < 
c||V 3 z||. Prom this, noting that I^z = I^IhZ , we conclude 

P{Uh){I ( 2h * -I ( 2h Z h) = h 2 p{u h )(l!$w) + h 3 p{u h ){T h ), 

and, using Galerkin orthogonality, 

p( u h)(4h Z ~ I 2 h z h) = h 2 p(u h )(I^w-I h w) + h 3 p(u h ){T h ). 

Assuming that the interpolation operators and 1^ behave like the H l -stable 
Clement operator, we have 

\I^w-I h w\\ K + \h} /2 \\I$w-I h w\\ dK < c/i||Vu/||s K . 

Collecting the above estimates, and using again (5.10), we find 

\p( u h){l 2 h z -I ( 2 h z h )\ < c(u, z)h 3 . (5.12) 

Finally, inserting the estimates (5.11) and (5.12) into (5.9), we conclude the desired 
estimate 

J(e) = E(u h ) + 0(h 3 ). (5.13) 

We emphasize that this estimate has been derived for a smooth dual solution z 
which excludes almost all interesting applications. Furthermore, the meshes are 
required to be uniform which conflicts with our ultimate goal of mesh adapta¬ 
tion. Therefore, a really meaningful analysis has to deal with the following three 
complications: 

• The estimates are to be proven on locally refined meshes with only limited 
uniformity properties. 

• The estimates are to be localized in order to allow for singularities in the 
dual solution. 

• Cases of non-smooth solutions u , either due to data irregularities or due to 
reentrant corners (the more severe case) should be considered. 
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with a cell patch K around K , if u is not regular (e.g., corner singularities) with 
constants independent of h . The first relation can be proven for quasi-uniform 
meshes, but is still open for general locally refined meshes including hanging nodes. 
We want to check these conditions by numerical experiments. 

a) Consider the Poisson model problem 

-A u = f in u\ 9 q = 0, 

on the domain Q = (— 1,1) 2 with f(x) = \' K 2 cos(^xi7r) cos(^X27t) . Apply the 
duality-based mesh refinement strategy for computing the two different functional 
values: 

J{u) = diu(0.5,0.5), J(u) = J d\u(l,X2) dx2, 

and monitor the behavior of D\u^x for increasingly refined meshes. 

b) Consider the Poisson model problem 

-A u = f in fi, u\ d n = 0, 

on the slit-domain fi = (—1, l) 2 \ {x € R 2 ,Xi =0,-1 < X 2 < 0} again with 
f(x) = cos(|xi7r) cos(|x27t) . Apply again the duality-based mesh refinement 
strategy to compute the same functional values as above, and monitor the behavior 
of the quantity 


Halloo := max {r^ 2 h^ l2 \\\d n u h }\\ K , }, 

^ € 11 h 

which reflects the expected singular behavior of u at the tip of the slit, where rx 
is the distance of the cell K from this point. Interpret the observed results. 



Chapter 6 

An Abstract Approach for 
Nonlinear Problems 


In this chapter, we will present a very general approach to a posteriori error estima¬ 
tion for the Galerkin approximation of nonlinear variational problems as developed 
in Becker and Rannacher [31]. The framework is kept on an abstract level in order 
to allow later for a unified application to rather different situations, such as non¬ 
linear PDEs, but also eigenvalue and optimization problems. We prepare for this 
by recalling the special situation of a linear Galerkin approximation of the form 


a(u,'0) = (6.1) 

a(u h ,iph) = F(iph) VVv» € Vfc, (6.2) 

with a linear functional F(-) as right-hand side. For a given (linear) output func¬ 
tional J(-) the associated continuous and discrete dual problems are 

a(ip , z) = J W>) Vty € V, (6.3) 

a(iph, z h) = J(^h) Viph.€V h . (6.4) 


Then, using Galerkin orthogonality, for the ‘primal’ error e u — u /, and the 

‘dual’ error e* := z — , there holds 

J(e) = a(e,z ) = a(e,e*) = a(u,e*) = F(e*), 

and, introducing the weighted primal and dual residuals, 

J(e) = F(z-ip h ) - a(u h ,z—i/> h ) =: p(u h )(z—i/) h ), %j) h € V h , (6.5) 

F(e*) = J((u-<p h ) - a(u-ip h ,z h ) -■ p*(z h )(u-<p h ), Ph&V h . (6.6) 

This trivially implies the ‘linear’ error representation 

J(e) = \p{u h ){z-ip h ) + \p f {z h )(u-ip h ), iph,i>h € V h . 


(6.7) 
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Below, we will see that this identity has a natural generalization to nonlinear 
problems, where the weighted primal and dual residuals are not identical. 


6.1 Galerkin approximation of nonlinear equations 

The following theory will be presented within an abstract functional analytic 
framework making as little assumptions as possible. Let A(u)(-) be a semilin- 
ear form and J(-) an output functional, not necessarily linear, defined on some 
function space V. The goal is the evaluation of J(u) from the solution of the 
variational problem 


A(u)(ip) = 0 VipeV. (6.8) 

The corresponding Galerkin approximation uses finite dimensional subspaces Vh C 
V to determine Uh G Vh by 

A(u h )(tl> h ) = 0 Viph G V h . (6.9) 

We assume the existence of directional derivatives of A and J up to order 
three denoted by A f (u)(<p, •), A"(u)(ip, p, •), A'"(u)(£, ip, p, •), and J'(u)((p ), 
J"(u)(ip, <p) , J'"(it)(£, ip, <p ), respectively, for increments G V. In these 

forms the dependence on the first argument in parentheses may be nonlinear while 
the dependence on all further arguments in the second set of parentheses is linear. 

Example 6.1. A typical example of a nonlinear problem of the form we are inter¬ 
ested in is the so-called vector Burgers equation 

—i/Au 4- u-Vu = /, in Q, u\qq = 0. 

for a vector function u G i/g(^) d • This is the natural generalization of the clas¬ 
sical 1-D Burgers equation to multiple dimensions. The corresponding variational 
formulation has the form (6.8) with the semilinear form 

A(u)('0) := i/(Vu, + (tt-Vu, ip) — (/, ip). 


Example 6.2. Another example of a nonlinear problem which will be considered 
in several exercises below is the diffusion-reaction equation 

-A u -u 3 = /, in fi, u\ d n = 0. 

for a scalar function u G Hq(Q) . This problem is interesting in the context of 
bifurcation theory. In this case the corresponding semilinear form is given by 


A(u)(ip) := (Vu, Vi/>) - (u 3 ,i>) - 
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For estimating the error J(u)—J(uh ), we employ the Euler-Lagrange method 
of constrained optimization. Introducing a ‘dual’ variable z e V (‘adjoint’ variable 
or ‘Lagrangian multiplier’), we define the Lagrangian functional 

C(u,z) := J(u) — A(u)(z ), 

and seek for stationary points {u, z} £ V x V of £(•, •), i.e., 

C'{u,z)(<p,i/>) = | ^a(u)(i/;) A | V{y>,V>}. (6.10) 

Clearly, the u-component of any such stationary point is a solution of the original 
problem (6.8). The Galerkin approximations {iih,Zh} € 14x14 are defined by the 
discrete Euler-Lagrange system 

C'(u h ,z h )( Vh ,4> h ) = { A,(Uh)(iPh,Zh)} }=0 VfaM, (6.11) 

where, again, the u/rcomponent of any stationary point is a solution of the discrete 
problem (6.9). Now, the goal is to estimate the error J(u)—J(uh) in terms of the 
residuals associated with this set of equations. We prepare for this by considering 
first the general situation of the Galerkin approximation of stationary points of 
functionals. 

Proposition 6.1. Let L(-) be a three-times differentiable functional defined on a 
(real or complex) vector space X which has a stationary point x G X, i.e., 

L f {x)(y)= 0 VyeX. (6.12) 

Suppose that on a finite dimensional subspace Xh C X, the corresponding Galerkin 
approximation 


L'(x h ){yh) =0 Vy/, e X h , (6.13) 

has a solution Xh € Xh • Then, there holds the error representation 

L(x) - L(x h ) = \L'(x h ){x-y h ) + TZ h , y h € X h , (6.14) 

with a remainder term IZh which is cubic in the error e x := x — Xh, 

1Zh •= \ [ L fff (xh+se x )(e x ,e x ,e x ) s(s— l)ds. 

Jo 

Proof. Since L f (x)(e x ) = 0, we can write 


L(x) - L(xh) = [ L'(x h +se x )(e x )ds 

Jo 


+ \L'{x h ){e x ) -\{L'{x h ){e x ) + L'(x)(e x )}. 
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The last term on the right is just the approximation of the integral term by the 
trapezoidal rule. For this, we have the well-known error representation 

f f(t)dt = §{/(0) + /(1)} + \ f f"(s)s(s—l) ds. 

Jo Jo 

Hence, we obtain 

L(x) - L(x h ) = ±L'(x h ){e x ) + \f L m {x h +se x )(e x , e*, e x ) s(s-l) ds. 

Jo 

Finally, observing that L f (xh)(yh) = 0 for all yh € Xh , we have that 

L f (x h )(e x ) = L f (x h )(x-y h ), y h G X h . 

This completes the proof. □ 

As an immediate consequence of Proposition 6.1, we obtain the following 
result for the Galerkin approximation of variational equations. 

Proposition 6.2. For any solution of equations (6.8) and (6.9), we have the error 
representation 

J(u) - J(u h ) = \p{u h ){z-tph ) + \p*(uh,z h )(u-iph ) + ^i 3) , (6.15) 

with arbitrary iphi'&h € 14, and the [primal ’ and ‘duaV residuals 

p(uh){-) ■■= -A(u h )(-), 

P*(uh,Zh){-) := J'{u h )(-)-A'(u h ){-,z h ). 

The remainder term is cubic in the ‘primal’ and ‘dual’ errors e u — Uh 

and e* := z — Zh , 

= h f {J"'(u h + se)(e,e,e) - A'"(u h +se){e,e,e,z h +se*) 

Jo 

— 3A"(uh+se)(e, e, e*)} s(5 — 1) ds. 

Proof. In order to apply Proposition 6.1, we define the space X := V xV and for 
arguments x = {u,z} G X the functional L(x) := C(u,z ). In this context the 
stationary points are denoted by x := {u, z} and Xh := {uh,Zh} with the error 
e x := x — Xh . Then, we have 

J(u) - J(u h ) = L(x) - A(u)(z) - L(x h ) 4- A h (u h )(z h ) = L(x) - L(x h ). 
Hence, the error representation of Proposition 6.1 gives us 

J(u) - J(u h ) = \L'(x h ){x-y h ) + 7 Z h , y h € X h , 
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with 

7 Z h := ^ f L'"(xh+se x )(e x ,e x ,e x ) s(s — l) ds. 

Jo 

By construction, we have for arbitrary yh = {(phi^Ph} € Xh : 

L'{x h )(x-y h ) = C' u (u h , z h )(u-ip h ) + C' z (u h , z h )(z—ip h ) 

= J'{u h ){u-tph) ~ A'{u h ){u-(ph,z h ) - A(u h )(z-iph) 

= p*(u h , Zh)(u-<Ph) 4- p(u h )(z-iph)- 


Notice that C(u, z) is linear in z. Consequently, the third derivative of L(-) 
consists of only three terms, namely, 


J"'(w/ l +se)(e,e,e) - A n, (uh+se)(e, e, e, Zh+se*) - 3A”(uh+se)(e, e, e*). 


This implies the asserted form of the remainder term 


(3) 



Remark 6.3. The derivation of the error representations (6.14) and (6.15) does 
not require the uniqueness of solutions; this is important, for example, for the 
application to eigenvalue problems. In cases with non-unique solutions, the a priori 
assumption xu —► x (h —► 0) makes the result meaningful as then the remainder 
term can be assumed to be small. 

Remark 6.4. The actual evaluation of the error identity requires guesses for the 
primal and dual solutions u and 2 which are usually computed from the Galerkin 
approximations Uh and zh by some post-processing as described in Section 4.1. 

Remark 6.5. The cubic remainder 7can usually be neglected. However, in 
parameter-dependent problems when approaching a bifurcation point, the deriva¬ 
tives of A(u)(-) and consequently 1Z^ may become large. In such a situation 
the abstract theory can still be applied but with special care. The extreme situa¬ 
tion will be seen in Chapter 7 in the context of eigenvalue problems where one is 
directly working in the singular case. 

In the linear case, we have seen that the primal and dual residual terms 
coincide, i.e. p(uh){z—'iph) = p*(zh){u—(fh) • This is no longer true in the nonlinear 
case, but the deviation from this property, i.e. the degree of nonlinearity of the 
problem, can be estimated as the following proposition shows. 


Proposition 6.6. With the notation from above, there holds 


p*{uh,z h ){u-ip h ) = p(u h )(z-tph ) + a p, 



for any <ph, iph € Vh , with 

A p:= / {A"(u h +se)(e,e,Zh+se*) - J"(iih + se)(e,e)} ds. 

Jo 
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Further, we have the simplified error representation 

J(u) - J(u h ) - p(u h )(z-ip h ) + 
for any iph €.Vh, with the quadratic remainder 

:= [ {A”(u h + se)(e,e,z) - J”(uh+se){e,e)} sds 

Jo 

Proof. We introduce the scalar function g(-) by 

g(s) := J'(u h + se)(e) - A f (u h + se)(e,z h +se*). 

By the definition of z and zh , there holds 

0(1) = «/'(w)(e) - A'(u)(e, z) = 0, 

0(0) = J'(u h )(e) - A'(u h )(e, z h ) = p*(u h , z h )(e), 

and 



g\s) = J”(u h +se){e,e) - A"(u h +se)(e,e,z h +se*) 


- A'(u h + se)(e,e*). 

Therefore, using Galerkin orthogonality, 


p*(u h , z h )(u-(p h ) = p*(u h , z h )(e) = g( 0) = p(0) - g(l) 


/Vw 

Jo 


ds 


/ {A (u h +se)(e,e, Zh+se*) - J (uh+se)(e,e)} ds 

Jo 

L 

+ / A'(uh+se)(e,e*)ds 
0 

= Ap + p(u h )(e*) = A p + p{u h ){z-\l) h ). 

this proves (6.16). In order to prove (6.17), we use integration by parts, obtaining 



n 


( 2 ) 

h 


/ {A"(uh + se)(e,e,z) - J"(uh+se)(e,e)} sds 

Jo 


= - {A'(u h +se)(e,z) - J'(u h +se)(e)} ds + A f (u)(e,z) - J'(u)(e), 

Jo 

where the last two terms vanish by definition of z. Consequently, employing again 
Galerkin orthogonality, we obtain 

= -p{u h )(z-jp h ) + j(u) - j(u h ), 

for arbitrary iph G Vh . This completes the proof. 

We note that the simplified error representation (6.17) could have been de¬ 
rived also from (6.15) using the relation (6.16). However, this involves lengthy 
calculations, so that we preferred to present a more direct argument. □ 
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In order to use the error representations (6.15) or (6.17) for practical mesh 
adaptation, we have to evaluate the primal and dual residual terms. As in the 
linear case, this requires approximation of the dual solution z and in the context 
of (6.15) additionally that of the primal solution u. This may be achieved again 
by post-processing of the Galerkin solutions Zh and Uh exploiting higher-order 
interpolation. Let the resulting approximations be denoted by Zh and Uh . Then, 
neglecting the remainder terms the approximate error representations take the 
form 


E(u h ,z h ) := \p(u h ){ih-z h ) + \p*(uh,Zh)(uh~Uh), 



and 


E(uh) := p(u h )(zh — z h ). (6.19) 

Remark 6.7. The identity (6.16) is useful as it offers the possibility of controlling 

the remainder 7in the simplified error representation (6.17). In fact, comparing 
the two error representations (6.15), (6.17) and using (6.16), we see that 

= -p(Uh)(z- 1 ph) + J(u)~ J{u h ) 

= -p{uh)(z-rph) + \p{uh){z-iph) + \p*{uh,z h )(u-<p h ) +n < £ ) 

= \p*{uh,z h ){u-p h ) - \p(u h ){z-il> h ) + 

= iA p + nf. 

Hence, we may try to control the linearization error by a posteriori checking the 
condition 


|Ap| « \p*(u h ,z h )(uh-u h ) - p(u h )(z h -z h )\ -C TOL. (6.20) 

where Uh~u and Zh ~ z are higher-order approximations and the cubic remain¬ 
der term 7is neglected. 

Remark 6.8. The possibility of improving on the linearization error , i.e. reducing 

the remainder term 7, by post-processing the Ritz approximation Uh has 
been studied in Vexler [134]. This costly process may be relevant in cases when 
the problem to be solved is close to bifurcation. 

Remark 6.9. A posteriori error estimates for the Galerkin finite element approxi¬ 
mation of nonlinear variational problems have also been derived using extensions 
of the classical ‘energy-norm-based’ approach; see, e.g., Verfiirth [131, 133]. These 
results rely on assumptions on the monotonicity of the underlying problem, i.e. 
the coercivity of certain derivative forms, or involve nonlinear stability constants 
which depend on the unknown solution. These estimates usually represent the 
‘worst case’ scenario and will only in special cases be of practical value. 
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6.2 A nested solution approach 

For solving the nonlinear problems by a Galerkin finite element method, we employ 
the following iterative scheme. Starting from a coarse initial mesh To , a hierarchy 
of refined meshes 

T 0 C Ti C • • • C T/ C • • • C T l , 

and corresponding finite element spaces VJ, l = 0, ...,L, with dimensions Ni is 
generated by the following nested solution process: 

1. Initialization: For l = 0, compute a solution uq £ Vo on the mesh To . 

2. Defect correction iteration: For l > 1, start with = w/_i £ VJ. For a 
computed iterate u\^ £ V\ evaluate the defect 

€ Vj, 

and solve the correction equation (Newton update) 

A'(u\ 3) )(v \ 3) , %pi) = , ipi) VV>i e V;, 

by Krylov-space or multigrid iterations using the hierarchy of previously 
constructed meshes {T/_i,..., To}. Update = u\^ with a 

step-length parameter , set j = j + 1, and repeat the iteration. This pro¬ 
cess is carried until a limit ui £ Vi , is reached with some required accuracy. 

3. Error estimation: Solve the (linearized) discrete dual problem 

z t eVi : A'(ui)((pi,zi) = J(<pi) VpieVt, 
and evaluate the a posteriori error estimate (6.18): 

J(e t ) « E(ui,zi). 

If \E(ui,zi)\ < TOL, or Ni > N max , then stop. Otherwise cell-wise mesh 
adaptation yields the new mesh T/+i. Then, set l := l + 1 and go to (2). 

Remark 6.10. The nonlinear iteration described above is oriented at the solution 
of stationary elliptic problems in which a global transfer of information is present. 
In the case of transport-dominated problems, particularly those with information 
transfer into one direction only such as in nonstationary problems, one would or¬ 
ganize the solution process differently, taking this transport direction into account. 

Remark 6.11. In the described Newton-like iteration the mesh adaptation is done 
for the limit solution on the current mesh in order to have a rigorous theoretical 
basis. However, it may be inefficient to carry the iteration on a coarser mesh to 
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the limit knowing that the discretization accuracy on this mesh is still insufficient. 
Hence, one would like to combine the estimation of the discretization error with 
that of the iteration error, both in accordance with the accuracy in the target 
quantity. Such a combined error estimator has been developed for a linear multigrid 
iteration in Becker et al. [27] and Becker [19], but it is an open questions for the 
nonlinear Newton iteration. 

Remark 6.12. The solution of the linear dual problem usually requires much less 
work compared to solving the nonlinear primal problem (6.8). In fact, in the con¬ 
text of the nested solution method described above the primal solution is obtained 
by a Newton iteration which requires the solution of several linear problems until 
convergence is reached. Then, solving the linear dual problem for the converged 
primal solution normally corresponds to about one additional Newton step. Fur¬ 
ther, this extra work is spent on optimized meshes adapted to the particular goal 
of the computation. Hence, particularly for nonlinear problems, the duality-based 
approach to adaptivity becomes relatively ‘cheap’. This is demonstrated by the 
examples from structural and fluid mechanics which will be presented in Chapter 
10 and Chapter 11. 

Remark 6.13. We remark on the following particular aspect in the evaluation of the 
dual residual p*(uh, Zh)(u—(ph) • The discrete dual solution Zh is determined by the 
variational equation (6.11) which is the Galerkin discretization of the correspond¬ 
ing continuous variational problem (6.10). The latter can usually be interpreted 
as the ‘weak’ form of a certain system of differential equations. The cell and edge 
residuals and of the dual solution zh are then taken with respect to the 
‘strong’ form of this dual problem. However, in deriving these residuals, we have 
to observe their consistency with the variational formulation (6.11) in order not to 
miss terms which may not be present on the continuous level, such as for example 
the residual terms occurring in solving the incompressible Navier-Stokes 

equations. We will comment on this point later on when a ‘naive’ derivation of the 
dual cell and edge residuals may lead to wrong results (see Remark 11.2). 


6.3 Exercises 

Exercise 6.1. Consider the Galerkin approximation of the stationary 1-dimensional 
Burgers equation 


-vu xx + uu x = /, in (0,1), u( 0) = 0 = u(l), 

by using piecewise linear finite elements. Determine the cubic and quadratic re¬ 
mainder terms in the a posteriori error representations for a linear output func¬ 
tional. 

Exercise 6.2. Use the second-order ‘nonlinear’ a posteriori error estimate to derive 
the L 2 -norm a posteriori error bound for the Galerkin finite element approximation 
of the (linear) Poisson problem. 
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Exercise 6.3 (Practical exercise). Consider the nonlinear diffusion-reaction equa¬ 
tion 

-A u - u 3 = /, in u\dn = 0, 

where SI is the domain defined in Exercise 3.3. The discretization is by the usual 
bilinear finite elements. For / = a = 1,..., 75, compute the solution’s mean value 
over the subdomain Q' = {x € X\ > .5, £2 > -5} C 

J\{u) := 4 / u(x) dx. 

Jet 1 

Monitor the size (relative to the estimated error) of the control term for the lin¬ 
earization error, 


A p := p*(u h ,z h ){I^u h -u h ) - p(u h )(l 2 h Zh-Zh), 

by approximating the weights using the biquadratic patch-wise interpolations 
l!$u h and I^h z h of the discrete primal and dual solutions Uh and Zh , respec¬ 
tively. Repeat the same calculations using the true biquadratic Ritz projection z 
of z instead. Explain the observed results. 



Chapter 7 

Eigenvalue Problems 


In the following, we will apply the abstract theory of the DWR method developed 
in Chapter 6 to error control in the approximation of eigenvalue problems. We 
mention some prototypical examples we are particularly interested in: 

• The symmetric eigenvalue problem of the Laplace operator: 

—A u = Xu. 

• The nonsymmetric eigenvalue problem of a convection-diffusion operator: 

—A u + b-Vu = Xu. 

• The stability eigenvalue problem governed by the linearized Navier-Stokes 
equations: 

—uAv + v-Vv + v-Vv + Vp = Av, V'V = 0, 

where v is some ‘base solution’ the stability of which is to be investigated. 

The results we will present for the approximation of these problems are taken from 
Heuveline and Rannacher [77, 78]. 

The above eigenvalue problems are all treated within an abstract setting, 
which will be laid out in the following. Let V be a (complex) Hilbert space. We 
seek {u, A} G V'xC satisfying 

Au — X Mu, 

with linear operators A and M. in V. The operator M. is assumed to be self- 
adjoint and positive semi-definite. It is introduced to allow for the situation of the 
stability eigenvalue problem for the Navier-Stokes equations, where the eigenvalue 
term only occurs in the first of the two equations (the momentum equation). We do 
not go deeper into the abstract functional analytic setting of eigenvalue problems 
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since we will concentrate below on concrete problems formulated in the standard 
function spaces. Again, we prefer the variational formulation 

a(u, t/>) = Vt/> £ V, (7-1) 

where a(*, •) := (A-, •) is the (generally nonsymmetric) sesquilinear form generated 
by A and ra(-,-) is the symmetric semi-definite sesquilinear form corresponding to 
M (usually the L 2 scalar product). Eigenfunctions are assumed to be normalized 
by m(u , u) = 1. For nonsymmetric A , one also considers the associated adjoint 
eigenvalue problem for the Hilbert-space adjoint A* of A , 

A*u* = \* Mu*. 

Observing that a*(u* 9 (p) = (A*u*,<p) = a(ip,u*) and m(u* 9 (p) = m((p 9 u*) 9 this 
has the variational form 


a((p,u*) = A*m(</?,u*) V</?€V. (7.2) 

Note that in the context of matrix eigenvalue problems the dual eigenvectors u* 
are also called left eigenvectors. From the definition, we see that primal and dual 
eigenvalues are related by A* = A, while the corresponding eigenvectors may 
differ. Usually, the dual eigenvectors are normalized by requiring 

m{u,u*) = 1. 

This normalization is possible only if u* is not m-orthogonal with respect to the 
whole eigenspace of A, which is equivalent to requiring that A has trivial de- 
fect , i.e., its algebraic eigenspace reduces to the geometric one (only trivial Jordan 
blocks in the matrix case); for a more detailed discussion of this issue see Heuve- 
line and Rannacher [77], and the literature cited therein. We will see that the 
simultaneous consideration of primal and dual eigenvalue problems is essential for 
rigorous a posteriori error estimation. 

The Galerkin approximation of (7.1) and (7.2) uses finite dimensional sub¬ 
spaces Vh C V 9 as described above, to determine A^}, {u* h , K) e v h xC, 
satisfying 

a(u h ,iph) = Km(uh,tph) € V h , m{u h ,u h ) = 1, (7.3) 

a(<Ph,K) = Km(<Ph,u* h ) Vi PhZV h , m(u h ,u* h ) = 1. (7.4) 

Our goal is to control the errors A —A/,, it. —«/,. and u* —u* h in the eigenvalues 
and eigenfunctions in terms of the residuals associated with these equations. 


7.1 A posteriori error analysis 

In order to derive a posteriori error estimates, we embed the present situation 
into the general framework of variational equations as considered in Chapter 6. 
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Consider the product spaces V := V xC and Vh := Vh xC C V and define for 
pairs U := {u, A} G V and = {xp, v} G V the semilinear form 

A(U)(^) := A m(xi, xp) — a(u, xp) +17 {m(u, u) — l}. 

With this notation, the above continuous eigenvalue problem and its discrete ana¬ 
logue can be written in the following compact form of semilinear variational equa¬ 
tions: 


A(U)(V)= 0 V^GV, (7.5) 

A{U h )(V h ) = 0 V* h g V fc , (7.6) 

where U := {it, A} and Uh := {w/i, A^} , respectively. In order to control the error 
in the approximation of the eigenvalues, we use the output functional 

J($) := <p). 

Since m(u, u) = 1 at the solution U = {w, A} , there holds 

J({7) = A m(u,u) = A, 

i.e., as desired this functional picks out the eigenvalue. We recall the Lagrange 
approach from Chapter 6, particularly the Lagrangian functional for arguments 
U — {u, A} and \I> = {xp, v} : 

C(U,V) = J(U)-A(U)(V) 

— A m(u,u) — A m(u,xp) + a(u, — z7 {m(u, u) — l}. 

The dual solution Z = {z,7r} G V and its Galerkin approximation Zh = 
{ 2 ^, 717 *} G are then determined by the equations 

j4'(17)($, Z) = J\U){$) G V, (7.7) 

^'(^)(^, Zh) = J'(U h ){$h) G V h . (7.8) 

The left and right hand sides of (7.7) read, for Z — {z, 7r} , U = {tx, A}, and 
$ = {ip, p} , as follows: 

A / (t/)(^>, Z) = Ara(<p,z) — a(<p, 2 )-b/xra(w, 2 :) + 27rRera((p,7x), 

«/'([/)($) = pm(u, u) + 2ARem(tp, u). 

Hence, the continuous dual problem takes the form 

Am(<p, z) — a(p, z) + p{m(u, z) — m(u, tx)} + 2{7r — A}Rera(<p, xx) = 0, 

for all $ = {ip, p} . We now show that the dual eigenpair U* = {xi*, A*} solves this 
equation, i.e., Z = {z, 7r} = U* is one solution. This is clear since for m(u,z) = 
m(u , it) = 1 and 7f = A , the equation reduces to 


a(<p, z) = 7T m(ip, z), 
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which is the defining equation for U *. Whether this is the only solution is not 
relevant for the present discussion. Analogously, the discrete adjoint problem (7.8) 
is solved by the discrete dual eigenpair Ufc = {w£, A£} determined by 

a(Vh,K) = K m (‘Ph,u*h) V<Ph€V h , m(u h ,u* h ) = 1. (7.9) 

With these preparations, we can formulate the following proposition: 

Proposition 7.1. With the primal and dual eigenvalue residuals 

p(uk, Afc)(-) := a(u h , •) - A h m(u h , ■), 

P*M, K)(-) :=a(;u* h )-\* h m(;u* h ), 
we have the error representation 

A-A* = \p{u h ,\ h ){u*-xp h ) + \p^ul,\* h ){u-<p h ) + n h , (7.10) 

for arbitrary t/^, iph £ Vh, with the cubic remainder term 

K h = |(A-A h )m(u-u h ,u*-u* h ). 

Proof. The assertion is an immediate consequence of Proposition 6.2 applied to 
the present situation. We have 

J(U)-J(U h ) = A {J'(U h )(U-* h )-A'(U h )(U-* h ,Z h )} 

+ ±{-A(u h )(z-* h )} + n h , 

for arbitrary $/, = {y/,.//}. 'I'/, = {iph,x} € V/,, with the cubic remainder term 
IZh which we evaluate below. Hence, using the above preparations, 

A-A h = p) m{uh,Uh) + 2X h Re m(u-iph,Uh) 

+ a(u-tp h , Zh) - A h m(u — (ph, z h ) - (A -p) m(u h , z h ) 

- 2nh Re m(u-<p h ,u h ) - ( X-p){m(u h ,u h )-l }} 

- i{A h m(u h ,z-iph) - a(uh,z—iph) + Or—x) {m(u h ,u h )-l}} 

+ Kh- 

We are free to choose p := A, x ~. and observing that m(u,u) = m(iih,Uh) = 
1, and Xh = 7th = A£ , we obtain 

A —Ah = X{a(u-tp h ,z h ) - X h m(u-<p h ,Zh)} 

+ \{a{uh,z-i)h)-X* h m{uh,z—ip h )} + Kb¬ 
it remains to evaluate the remainder term. Setting E := {u — Uh,X — A*} and 
E * := { u* — u* h . A* — \* h }, the general remainder term from Proposition 6.2 is 

K h = A [ { J" , (U h +sE)(E,E,E)-A , "(U h +sE)(E,E,E,Z h +sE *) 

Jo 

- 3A"{U h +sE){E, E, E*)}s(s-1) ds. 
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In the present case, by a simple calculation, we have 

J ff '{Uh+sE)(E, E, E) = 6(A-A h )m{u-u h ,u-u h ), 
A'"{U h +sE){E,E,E,Z h +sE *) = 0, 

-3A"(U h + sE)(E, E, E*) = -6(A-A h )m(u-u h ,u*-u* h ) 

- 6(\*-\* h )m(u-u h ,u-uh)- 

Consequently, noting that A —A^ = A* — A£ , 

IZh = — 3 / (A — \h) rrt(u—Uh, u* — Uh)s(s — 1) ds 

Jo 

= |(A-A h )m(u-u h ,u*-u* h ), 

which completes the proof. □ 

Remark 7.2. We add the following remarks concerning Proposition 7.1: 

• The error representation (7.10) holds true without any assumption on the 
multiplicity of the eigenvalue A or its defect. However, such a restriction 
will become necessary below in dealing with the error of the eigenfunctions. 

• In the error representation (7.10) only terms involving the computed primal 
and dual eigenpairs occur and no additional outer dual problem needs to be 
solved. We will see a similar situation in the context of optimization problems 
discussed in Chapter 8 where the underlying mechanism will become clear. 

• In the nonsymmetric case the simultaneous consideration of primal and dual 
eigenvalue problems is essential within an optimal multigrid iteration anyway 
(see Heuveline and Bertsch [76]). Then, the computation of u* for the error 
estimator does therefore not introduce extra work. 

• In Proposition 7.1, we have assumed that the governing operator A re¬ 
mains unchanged under discretization, i.e. all coefficient are frozen. Below, 
in considering the stability eigenvalue problem, we will additionally allow for 
approximation of the operator A(u) depending on coefficients u that will 
also be subject to approximation. 

Practical evaluation of the error representation 

Next, we want to determine the explicit form of the residuals in Proposition 7.1. 
To this end, we need to be more specific about the particular structure of the 
eigenvalue problem considered. Here, we restrict ourselves to a simple model sit¬ 
uation which, nevertheless, is prototypical for the problems we are interested in. 
On a polygonal or polyhedral domain consider the eigenvalue problem of 

a second-order elliptic differential operator A such as, for example, 

Av —Av -f b-Vv = Xv in Q, vj^q = 0, 
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with a smooth (or even constant) transport coefficient b. In this case, we have 
M := id. Further, let this eigenvalue problem be approximated by the Galerkin 
method using piecewise linear or d -linear finite elements on meshes = {K} , as 
described in Chapter 3. Within this setting, we can proceed analogously as before, 
obtaining 

p(u h ,\ h ){-) = a(u h ,-) ~ X h m(u h ,') 

= {^ Uh ~ ^hMlLh ,-)k — {dnUh,-)dI<} 

KeT h 

= ^2 {('AUh-XhMUh, ')k + \{[®n u h], ')dK}, 

KeT h 

P*( u h x *h)(-) = a (‘, z h ) — X* h m(-, z h ) 

= ^2 {(')'A* z h~^hM z h)i< ~ {',d^Zh)dK} 

KeT h 

= ^ {('^* Z h~K^ z h)K + !(*,[&n z h])dx}- 
KeTh 

Hence, using again the notation of ‘equation’ and ‘jump residuals’ Rh , 

and r£, respectively, analogously as introduced in Section 3.1, the residual admits 

the estimate 

\p(u h ,\ h )(u*-tP h )+p*(ul\* k )(u- Vh )\ < Y, {pkUk+PkVk}, (7.11) 

KeTh 

with the cell residuals p/c, p* K and and weights uk, defined by 

Pk := (\\Rh\\ 2 K + Ik/iHl/c) ^ , 

Pk (\\Rh\\ 2 K + ^ ll r hlli/f) ^ > 

WK := {\\u-<Ph\\ 2 K + 5*k 2 ||w-^/.HIa-) 1/2 , 
w k := (ll u * — i’hW'ji + ||it* — VVilliif) ^ ■ 

As a consequence of the above discussion, we obtain the following result: 
Proposition 7.3. Within the above setting , assuming that 

m(u-u h ,u*-u* h ) | < 1, (7.12) 

we have the ‘ weighted’ a posteriori error estimate 

|A-A fc | < <:= {pk»k + PkUk}, (7-13) 

KeTh 
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and the ‘energy-norm-err or-type’ estimate 

|A —A^| < r)^ := c\ h 2 i<{p 2 K + P*k }» (7-14) 

KeT h 

with a constant c\ growing linearly with |A|. 

Proof. Using the estimate (7.11) in the error representation (7.10) gives us 

|A — Xh\ < \ {pk^k + Pkljk} + [Rhl- 

KeT h 

Since, in virtue of assumption (7.12), 

\T^h\ = \\{\-\h)m(u-u h ,u*-u* h )\ < \ |A-Afc|, 

the asserted estimate (7.13) follows. To prove (7.14), we choose iph := IhU* and 
(fh := IhU in an d 00 k > with the modified nodal interpolation operator Ih 
introduced in Section 3.2. Writing 

u*-I h u * = (u*-u* h ) - I h {u*-u * h ), u-~I h u = (u-Uh) - h(u-Uh), 
in the weights u>* K and u>k, we obtain by the interpolation estimate (3.10) that 

h^ 2 {wf + w 2 k } < c^{||V(u*-0|| 2 + ||V(u-u h )|| 2 }. 

K<=T h 

This gives us 

|A —Afc| < ic/( ^ /i^{ / 9 2 f +p^}) 1/2 (||V(n*-0|| 2 + ||V(n-tx h )|| 2 ) 1/2 . 

KeT h 

Now, the proof would be completed by showing that the energy-norm errors of 
the eigenfunctions are proportionally bounded by the eigenvalue error. This is 
actually the case but requires lengthy calculations employing duality arguments 
for the eigenfunction errors. These details are omitted and we refer instead to 
Heuveline and Rannacher [77]. □ 

Remark 7.4. The analogue of the a posteriori error estimator 77 ^ for symmetric 
eigenvalue problems has been given by Nystedt [108]. There, the symmetry of the 
problem is extensively used. Furthermore, // 2 -regularity of the eigenfunctions is 
required which excludes domains with reentrant corners, the most interesting case 

for an adaptive approach. Notice that our derivation of 77 ^ does not need these 
assumptions. 
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Remark 7.5. An alternative version of the eigenvalue-error estimator 77 ^ has been 
given by Larson [101], also for the symmetric case and assuming H 2 -regularity of 
the eigenfunctions, namely, 


+ • (7-15) 

KeT h 

Both error estimators, 77 ^ and 77 ^ , are asymptotically equivalent for regular 
situations. However, the required H 2 regularity for the eigenfunctions renders the 

estimator 77 ^ useless for typical situations in which mesh adaptivity is needed. 

Remark 7.6. Neglecting the presence of the dual eigenvalue problem, we may try 
to control the error in approximating the eigenvalue using the primal residual part 
alone, i.e. using the ‘reduced’ error estimator 

W T \ d ^kPk- 

KeT h 

Below, we will compare the performance of the above eigenvalue-error estimators 
77 ^ , 77 ^ , 77 ^ , and 77^ ed for a simple model situation. 

Numerical examples 

We consider the convection-diffusion model eigenvalue problem 

—Av + b-Vv = Xv in Q, v^n = 0, 

on the slit-domain Q = (—1, l)x(—1,3) \ {x € R 2 , X\ =0,-1 < X 2 < 0} , with the 
transport vector b = (0, b y ) T ; for a sketch of this configuration, see Figure 7.1. In 
the computations on this test problem the mesh refinement is organized according 
to the ‘fixed-rate’ strategy with refinement rate X = 0.2, see Section 4.2. 

Test 1: Symmetric case 

At first, we consider the symmetric eigenvalue problem, i.e. b = 0. The eigen¬ 
function corresponding to the smallest eigenvalue and the computational mesh 
generated on the basis of the weighted error estimator 77 ^ are shown in Figure 
7.1. Figure 7.2 shows the mesh efficiencies achieved on the basis of the different er¬ 
ror estimators 77 ^ , 77 ^ and 77 ^ introduced above compared with that of uniform 

mesh refinement. (In this case 77^ ed and 77 ^ are equivalent.) We see that in the 
symmetric case all estimators show almost equally good performance compared to 
that of uniform refinement. This is due to the dominance of the error caused by 
the slit singularity which is well captured by all residual-based error estimators. In 
fact, it is known that in the symmetric case, the eigenvalue error is proportional 
to the square of the energy-norm error. 
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Figure 7.1: Configuration for b = 0 (left), adapted mesh with 12,000 cells (mid 
die), eigenfunction (right); from Heuveline and Rannacher [77]. 
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Figure 7.2: Mesh efficiency achieved using the different error estimators, rj\ 

( 2 ) 

(symbol o), rj\ (symbol U), and rfi) (symbol *), compared against uniform 

refinement (symbol A). The curves for and rf£ lie above each other; from 
Heuveline and Rannacher [77]. 
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Test 2: Nonsy] 



metric case 


Next, we consider the nonsymmetric version of the test eigenvalue problem with 
vertical transport, b y = 3 . In this case, due to the Dirichlet boundary conditions, 
the primal eigenfunction has a steeper gradient at the top boundary, while the dual 
eigenfunction has one at the bottom boundary; see Figure 7.3. The latter boundary 
layer strongly interferes with the slit singularity. Figure 7.4 shows adapted meshes 

obtained using the ‘energy-type’ error estimator , the ‘reduced’ error estimator 

rfff d , and the ‘weighted’ error estimator rff . The superiority of the latter one is 
clearly seen in Figure 7.5. 




Figure 7.3: Primal (left) and dual eigenfunction (right); from Heuveline and Ran- 
nacher [77]. 




Figure 7.4: Adapted meshes with 10,000 cells on the basis of the error estimators 
G e ft)’ V r \ ed (middle), (right); from Heuveline and Rannacher [77]. 
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Figure 7.5: Mesh efficiencies achieved on the basis of the different error estimators: 

rj ^ (symbol o), rf[ ed (symbol X), and rff (symbol *), compared against uniform 
refinement (symbol A); from Heuveline and Rannacher [77]. 


7.2 Error control for functionals of eigenfunctions 

We now turn to the question of how to estimate the error with respect to out¬ 
put functionals of the eigenfunctions. This question seems non-trivial since at an 
eigenvalue the corresponding adjoint operator is naturally singular, so that it is 
not clear what the proper definition of the dual problem has to be in this case. 

Let j(-) : V —> C be an output functional, for simplicity assumed to be linear, 
with respect to which the error in the eigenfunctions is to be controlled. For ease 
of presentation, we also assume that the eigenvalue A considered is simple with a 
normalized eigenfunction u. In order to apply our abstract theory, we define the 
functional 

Then, the associated abstract dual problem for Z = {z, 7r} 6 V, 

A'(U)($, Z ) = J'(E/)($) e V, 


takes the explicit form 


Am((^, z) — a(<p, z) + pm(u, z) — 27r Re ra(<p, u) = 


(7.16) 
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for all $ = {ip, pj G V . This is equivalent to 

A m(ip, z) — a(ip, z) = j (ip) + 2 7f Re m(<p, u) G V, 

and m(u,z) = 0. Since A is assumed to be a simple eigenvalue, by virtue of the 
Fredholm alternative, this equation can be solved if and only if its right-hand side 
vanishes on the eigenvector u , that is 

j(u) — 2 7r Rem(u, u) = 0 <=> 7 t~^j(u). 

Consequently, for the given eigenpair {w, A} the so reduced dual problem 

a(p, z) — \m((p, z) = j(ip) — j(u) Re m(ip, u) V<p G V, (7.17) 

has a solution 2 G V, which is uniquely determined in view of the property 
m(u,z) = 0. By an analogous argument, the reduced discrete dual problem is 
seen to be 


a(<Ph, z h ) - \ h m(<p h , z h ) = j(<p h ) - j(u h ) Re m((p h ,u h ) V(p h G 14, (7.18) 

where {uh,\h} is the eigenpair of the approximate eigenvalue problem, and 
ith •= 5 j{ u h) • This problem also has a solution Zh G Vh , uniquely determined by 
m(uh 1 Zh) = 0. We see that the well-posedness of these dual problems is guaran¬ 
teed by filtering out the eigenspaces span{w} and span{u/ l }, respectively. With 
all this, we can state the following proposition. 

Proposition 7.7. Let {uh,\h} be a computed eigenpair approximating {w,A}. 
Then, for the given functional j(-) : V —► C and the associated solution z G V of 
the dual problem (7.17), we have the error representation 


j(u u h ) = p(u h , X^iz-iph) + (\-\h)m(u-u h ,z) 

+ \j(u)m(u-u h ,u-u h ), 



for arbitrary iph G Vh . 

Proof. First, we recall the definitions of the primal (eigenvalue) residual 

p(uf 1 , A h )(-) = a(u h , •) - A h m(u h , •), 
and the dual residual associated to the dual problem (7.17), 

p*(u h ,z h )(-) := a(-,z h ) - \hm(-,z h ) + j{-) - j(u h )R£m(-,u h ). 

This implies with the results of Chapter 6 that 

j{u-u h ) = \ {a{u h ,z-iph) ~ ^hm(u h ,z-iph)} 

+ \ {a(u-(p h , z h ) - A h m{u-tp h , *h) 

-j(u h ) Re m(u-(p h ,u h ) + j{u-(fi h )} +Kh, 
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and consequently, taking ph = Uh , 

j(u-u h ) = a{u h ,z-iph) - ^hm(u h ,z-tph) + a(u-u h ,z h ) 

- A h m{u-uh, z h ) - j(uh) Re m{u-u h ,uh) + 2H h . 

To identify the remainder IZh , we note that 

J"'(U h +sE;E,E,E) = 0, 
A"\U h +sE-,Z h +sE*){E,E,E) = 0, 

-3A"(Uh+sE-,E,E,E*) = -6(A-A h )m(u-Uh, z-z h ) 

- 6(Tr-n h )m(u-u h ,u-Uh), 


which yields 

IZh - \{\-Xh)'m(u-Uh,z-Zh)-\-\{Tt-Tth)m{u-Uh,u-u h ). 

We recall that n = \j(u ) and nh = \j{uh) and obtain 

lZ h = i(A-A h )m(u-u h ,z-z h ) + \j(u-uh)m(u-uh,u-u h ). 

From this, we infer as an intermediate result: 

j(u-u h ) = a(u k ,z-ip h ) - \ h m(u h ,z-ip h ) 

+ a(u-u h , z h ) - A h m(u-u h , z h ) 

- j(u h ) Re m{u-u h ,u h ) + \ j{u-uh)m{u-uh,u-uh) 

+ (A-A h )m(u-Uh,z-z h ). 

Next, by definition and since m{uh,Zh) = 0, we have 

a(u-u h ,z h ) - \ h m(u-u h ,z h ) = (A-A h)m{u,z h ) = (A-A h )m{u-u h ,z h ). 

Further, noting that m(u, u) — m(uh,Uh) = 1, there holds 

m(u — Uh,u — Uh) = m(u,u) + m(^,%) - 2Re m(u,%) 

— m(u , u) — m{uh , Uh) — 2 Re m(u — Uh,Uh) 

= — 2 Re m(u — Uh, Uh)- 

Then, combining the last three relations gives us 

j(u~u h ) = a(u h ,z-ip h ) - \ h m(u h ,z-'ip h ) + (\-\ h )m(u-u h ,z h ) 

+ (it) 771 ( 14 - 1 ^, 11 - 1 ^) + (A-A/i)m(u-w/i,2J-2J/i), 

which completes the proof. □ 


Remark 7.8. The proposition requires A to be simple. In the case of geometric 
multiplicity p > 1, we have to simultaneously consider a whole basis , i — 
1, ...,p} of the eigenspace kern(*A — XI) in setting up the dual problem. The case 
of higher algebraic multiplicity can also be handled but is much more involved. 
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Numerical example 

In order to illustrate the foregoing result, we consider the model problem from 
above. The goal is to evaluate the derivative value j(u) := d\u(a) of the ‘first’ 
eigenfunction at a — (0.5,2.5) T . 



Figure 7.6: Mesh efficiency of rff (symbol *) compared to r)u ^ (symbol O) and 
uniform refinement (symbol A); from Heuveline and Rannacher [77]. 



Figure 7.7: Adapted mesh with about 12,000 cells for the approximation of d\u(a) 
obtained by rff (left), and the corresponding dual solution (middle and right); 
from Heuveline and Rannacher [77]. 
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7.3 The stability eigenvalue problem 

Eigenvalue problems play an important role in the analysis of the stability of 
solutions of nonlinear differential equations. Classical areas of applications are 
structural mechanics and fluid mechanics, in the latter context referred to as hy¬ 
drodynamic stability. We will consider this kind of problem again within an ab¬ 
stract setting. Let a(-)(*) be a given semilinear form determining the base solution 
ueV by 


a{u){xp) = 0 V-0GV, (7.20) 

and its Galerkin approximation Uh G Vh C V by 

a{uh)Wh)= 0 W'lph G V h . (7.21) 

Note that u is a stationary solution of the system under consideration. We want to 
determine whether this base solution is dynamically stable , i.e., whether any non¬ 
stationary solution trajectory {u(t) G V,t> 0 }, starting from a small perturbation 
fx(0) = u 4- w° , and satisfying the evolution equation associated to (7.20), 

m{d t u, ip) 4- a{u){iP) = 0 Wip G V, (7.22) 


stays bounded or even decays back to u, as t —► oo. Here, the decay is expressed 
in the semi-norm m(-,-) 1 / 2 , which in most practical cases is actually the usual 
L? norm. However, in order to prepare for the application of this approach also to 
the incompressible Navier-Stokes equations in Chapter 11, we allow ra(-, ♦) to be 
slightly more general. 

In the following, we outline the basics of the so-called linearized stability 
theory. The perturbation w(t) := u(t) — u satisfies the nonlinear perturbation 
equation 


0 = m(dtw, ip) + a(u+w) (ip) — a(u) (xp) 



m(d t w, ip) + a f (u)(w, xp) 4- / a''(u+sw)(w, w, ip)s ds, 

o 


for t > 0. Taking \p = w , we obtain 

+ a'(u)(iv, w ) 



a"(u-\-sw)(w, w, w)sds. 


o 


Assuming now that the perturbation w(t) is initially small such the the cubic term 
on the right can be neglected, the initial decay or growth of w(t) is determined by 
the spectral properties of the sesquilinear form a'(u)(-,-) obtained by linearizing 
a(fi)(-) at the base solution. Particularly, if the nonsymmetric stability eigenvalue 
problem 


a r {u){u^\p) — A m(u,\p) Vip G V , 


(7.23) 



96 


Chapter 7. Eigenvalue Problems 


has an eigenvalue X crit with negative real part, then any perturbation having 
initially a nontrivial component in the direction of an eigenfunction of X crti will 
grow initially and may eventually blow up. 

We want to solve this stability eigenvalue problem numerically by a Galerkin 
approximation which reads as follow: 

a'(uh)(uh,tph) = ^hm(u h ,iph) e V/>. (7.24) 

Our goal is to estimate the error in the critical eigenvalue, \ crtt — A£ r,t , i.e. the 
‘first’ eigenvalue with possibly negative real part, in terms of the residuals corre¬ 
sponding to primal and dual equations. Note that this includes both the error due 
to approximating the base solution u , as well as of the eigenvalue. To this end, 
we embed this situation into the general framework of variational equations laid 
out above. We introduce the spaces 

V := V x V x C, V h := V h x V h x C, 

A 

consisting of elements U := {u,u, A}, = {^,^, v} and Uh := {&/i, W/i, A/*}, 

A 

\f/ h = {i/j h , h , i/}, respectively. Using the semilinear form 

A(U)(^) := a(u)(^) 4- Am(u,^) — a f (u)(u,xp) + v{m(u, u) — l}, 
the continuous and discrete problems can be written in compact form as follows: 

A(U)(V) = 0 V^eV, (7.25) 

A(U h )(9 h ) = 0 Vyf'h G V fc , (7.26) 

For controling the eigenvalue error, we work again with the functional 

J($) := nm(<p,<p), $ = 

such that J(U) — X m(u, u) — X . Then, the corresponding continuous and discrete 
dual solutions Z = {i,z,7r} G V and Zh = {zh,Zh,ith} € are determined by 
the problems 

A'([/)($, Z) = J'(U)($) G V, (7.27) 

A'{U h ){* h , Z h ) = J'(t4)(*fc) G V fc . (7.28) 

A detailed computation shows that these equations are solved by Z = U* := 
{z,u*,A} and Z^ = U£ := {zh,u^,Xh}, respectively, where and are the 
(normalized) dual eigenfunctions as before, and the dual base solutions z and Zh 
are determined by the dual problems 

a!(u)(ip,z) = —a f \u)(tp,u,u*) V<p G V, 

a'(u h )(<p h ,z h ) = -a"(u h )(<p h ,u h ,u* h ) V(p € V. 


(7.29) 

(7.30) 
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The corresponding residuals of the approximate base solutions are 

■= -a(uh)(-), 

P*(uh,Zh)(-) ■= a"(u h )(-,u h ,u* h ) - a'(u h )(-,z h ), 
and those of the eigenpair approximations, 

p{u h , U h , A/j)(-) := a'(u h )(u h , ■) - X h m(u h , ■), 
p*(u h ,u* h ,\* h )(-) := a'(u h )(-,u* h ) - A * h m(-,u* h ). 

We collect the foregoing findings in the following proposition: 

Proposition 7.9. With the above residuals, we have the error representation 


A —A = 



\p(^h){z-li>h ) + , Zh)(u~<Ph) 

\p{u h ,u h , Xh){u*-iph) + \p*{uh,u* h ,X* h ){u-iph) 


+ Ti-hi 


(7.31) 


A 

for arbitrary 'ifth, ifth, <Ph, <Ph £ 14. The remainder 7 Zh is cubic in the errors 


e u 


u 


— Uh , e z := z — Zh , and e x := A —A^ , e u : 


u-u h , e 




u*-u* h : 




\e x m(e u ,e u *) 4- ^ 




//// 


(i4 4- se u )(e u , e u , e u ,Uh 4- se u , 4- se n *) 


- a"'(uh + 5c u )(c ,i , e tt , e u , + se*) 

4- 3 a”'(uh 4- se u )(e u , e u , e u , 4- se u *) 

+ 3a f "(uh + se u )(e u , e u , + se u , e“*) 

+ 6a"({4 + se“)(eVVn 

— 3 a ff (uh 4- se u )(e u , e u , e u *)}s(s - 


1 ) ds. 


Application to a model case 

We want to apply this abstract result to the concrete situation of the nonsymmetric 
model problem from Example 6.1 (vector Burgers equation). Here, the stability 
eigenvalue problem has the form 

-Au4u-Vu4u*Vw = Xu, in fl, u\dn = 0. (7.32) 

Its finite element approximation computes the discrete base solution Uh by solving 

(Vfifc, VVVi) + ( u h -Vu h , VVi) = (/, il>h) Wvi € 14, (7.33) 

and then determines the corresponding eigenvalue A h from the eigenvalue problem 
{Vu h , Viph) 4- (u h -Vu h +u h -Vu h ,ilj h ) = A h (uh^h) 4- (7.34) 
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The corresponding discrete dual problem determining Zh £ Vh reads 

(Vph, Vz h ) + (<p h ,Vu h z h -Viu h <8> z h )) 

= (ip h ,Vu h u* h -'V-(u h ®u* h )) V<Ph.€Vh, 

and that determining the discrete dual eigenpair {u£, A£} , 


(V<^, V<) + (<p h ,\7u h u* h -\7iu h ® <)) 
Then, the associated cell residuals are 


Ki<ph,<) v<phev h 


R 


h\K 


R*h\K 


R 


h\K 


D* 

n h\K 


f + Au h - Uh-Viih, 


A z h + Vu h u* h - V-(ufc <g> - Vu/jZ/, + V-(«h <g> Zfe), 

-Aufc + Wh-Vuh + u h -Vu h - X h u h , 


A tit 


V-(tih ® O + ViifcUfc - , 


(7.35) 


(7.36) 


while the associated edge residuals r h \ p, r^ r and r h j r , r^ r have the same form 
as above in the case of the simple Poisson equation. For this situation Proposition 
7.9 yields the following result: 

Proposition 7.10. Using the notation from above , and assuming again that 

m{u—uh,u*—u* h )\ < 1, 


we have the a posterior error estimate 

|A-Ah| < rjf := ^2 {pk&k + Pk&k + Pk^k + PkVk} + Hh- (7.37) 

Ker h 


The cell residuals and weights are defined by 

Pk := 

(M 2 * + ^ 1/2 ^llS*) 1/a , 

^ * 
lu k 


** 

• • 

II 


lj k := 

(P“ 0h dx) ! > 

Pk •= 

(ll^lli: + h- K 1/,2 ||r^||| K ) ^ , 

♦ 


♦ 

Pk : ~ 

(II^IIk + h K t || r* h q K ) ^ , 

wjc := 

( U — ifhW^K + ^K 2 U — VhWdK) ^ • 


for arbitrary iph, 0h, <Ph £ Vh, cmd the remainder IZh is cubic in the errors 

I A 7 AO I Of A U. A Of 

e u := u — Uh and e := u —u h : 

n h = - (V(e“ ® e“ - ±e“<g>e“), e u *). 
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An error estimator as derived in Proposition 7.10 will be applied below in Chap¬ 
ter 11 to the approximation of the stability eigenvalue problem of the Navier-Stokes 
equations. There, we will demonstrate the interplay of the different components 
in the error estimator . 

7.4 Exercises 

Exercise 7.1. For the Galerkin approximation of a symmetric eigenvalue problem 

a(u, ip) = A (u, ip) G V, 

with a symmetric bilinear form a(-, •) on a (real) Hilbert space V , the general a 
posteriori error representation reduces to the form 

A-A h = p(u h ,X h )(u-'ip h ) + |(A-A h )||w—«h|| 2 , € V h . 

Derive this identity by direct algebraic manipulation. 

Exercise 7.2. Consider the model convection-diffusion eigenvalue problem 

—Av 4- b Vv = Xv in vj^ = 0. 

Under the assumption that the eigenfunctions have H 2 - regularity, and that 

\(u-u h ,u* -u* h )\ < 1, 
prove the a posteriori error bound 

|A-Ah| < c x(^ ^2 + » 

Ker h 

with a constant ca = 0(|A|). 

Exercise 7.3. Consider the d-dimensional Burgers equation 

—vAu 4- u-Vu = / in 11, u^q = 0, 

for a vector function u : — > R d . Suppose that u G V := is a solution. 

Show that if all eigenvalues of the symmetric stability eigenvalue problem 

—vAw 4- ^{Vu4- Vu T — V ul}w = A w 

are positive, then u is dynamically stable. This criterion for stability is much 
stronger than that of linearized stability since it allows perturbations of any size 
and also applies to nonstationary solutions. Therefore, it cannot be expected to 
be satisfied in many practically interesting situations. 
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Exercise 7.4 (Practical exercise). Consider the nonlinear boundary value problem 
of Exercise 5.3, 

-A u - u 3 = /, in fi, u\ d Q = 0, 

where fl is the domain defined in Exercise 3.3. The discretization is by the usual 
bilinear finite elements. For / = a = 1,...,75, compute the solution on increas¬ 
ingly refined meshes and investigate its dynamic stability using the criterion of 
linearized stability. Monitor the behavior of the eigenvalue error estimator and, in 
particular, its components measuring the approximation of the base solution and 
that of the eigenvalue. Explain the observed results. 
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